DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 



Approved for public release; distribution is unlimited 



Design of a GPS Aided 
Guidance, Navigation, and Control System 
for Trajectory Control of an Air Vehicle 



by 



Eric N. ^lallberg 
Lieutenant, United States Navy 
B.S. University of Pennsylvania, 1984 



Submitted in partial ful fillm ent of the 
requirements for the degree of 

MASTER OF SCIENCE IN AERONAUTICAL ENGINEERING 



from the 



NAVAL POSTGRADUATE SCHOOL 



March , 1994 



Department of Aeronautics and Astronautics 



UNCLASSIFIED 
iliftITV dLA«IFI<:ATI6N 6F THIS P-AdE 



REPORT DOCUMENTATION PAGE 



REPORT SECURITY CLASSIFICATION 

UNCLASSIFIED 


lb. RESTRICTIVE MARKINGS 


SECURITY CLASSIFICATION AUTHORITY 


3. DISTRIBUTION/AVAILABILITY OF REPORT 
Approved for public release; 
distribution is unlimited. 


DECLASSIFICATION/DOWNGRADING SCHEDULE 


PERFORMING ORGANIZATION REPORT NUMBER(S) 


5. MONITORING ORGANIZATION REPORT NUMBER(S) 


NAME OF PERFORMING ORGANIZATION 
Naval Postgraduate School 


6b. OFFICE SYMBOL 
(if applicable) 

AA 


7a. NAME OF MONITORING ORGANIZATION 
Naval Postgraduate School 


ADDRESS (City, State , and ZIP Code) 

Monterey, CA 93943 


7b. ADDRESS (City, State, and ZIP Code) 

Monterey, CA 93943 


i. NAME OF FUNDING/SPONSORING 
ORGANIZATION 


8b. OFFICE SYMBOL 
(if applicable) 


9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 


ADDRESS (City, State, and ZIP Code) 


10. SOURCE OF FUNDING NUMBERS 


PROGRAM 
ELEMENT NO. 


PROJECT 

NO. 


TASK 

NO. 


WORK UNI ! 
ACCESSION* 



. TITLE (Include Security Classification) 



Design of a GPS Aided Guidance, Navigation, and Control System for Trajectory Control of an Air Vehicle 
^PERSONAL AUTHOR(S) 



LT. Eric N. Hallberg 



la. TYPE OF REPORT 


13b. TIME COVERED 




14. DATE OF REPORT (Year, Month, Day/ 


15. PAGE COUNT 


Master’s Thesis 


FROM 01/92 JO 


03/94 . 


March 1994 


119 



’* SUPPLEMENTARY notation views expressed in this thesis are those of the author and do not reflect the 

EAdal policy or position of the Department of Defense or the United States Government. 



\ COSATI CODES 


18. SUBJECT TERMS (Continue on reverse if necessary and identify by block number) 


FIELD 


GROUP 


SUB-GROUP 


Trajectory Control, GNC, Linear Quadratic Regulator, LQR, Kalman Filt 








Nonlinear Simulation, P- Implementation, Bluebird 









ABSTRACT (Continue onreverse if necessary and identify by block number) 



The advent of GPS has afforded the aerospace controls engineer a powerful, new means of controlling air vehic 
bis work explores a new method of designing and implementing controllers and guidance systems for autonomous con 
air vehicles utilizing a GPS integrated guidance, navigation and control system. This is a subject of consider* 
terest when realizing controllers to track reference trajectories given in an inertial reference frame. The desi 
iplementation, and dynamic simulation of a precise tracking trajectory controller for an Unmanned Air Vehicle (Uj 
presented. This design provides a natural conversion of commands and other measured outputs (such as GPS sign 
)m an inertial reference frame to a body-fixed reference frame. This achieves automatic recruiting of the actual 
iiile preserving the properties of the original design (linearization principle). 



i. DISTRIBUTION/AVAILABILITY OF ABSTRACT 
3 UNCLASSIFIED/UNLIMITED □ SAME AS RPT. □ DTIC USERS 


21. ABSTRACT SECURITY CLASSIFICA 

UNCLASS1 


TION 

[FEED 


a. NAME OF RESPONSIBLE INDIVIDUAL 

Isaac I. Kaminer 


22b. TELEPHONE (Include Area Code ) 

(408) 656-2972 


22c. OFFICE SYMBOL 

AA/KA 



All other editions are obsolete 



UNCLASSIFIED 



ABSTRACT 



The advent of GPS has afforded the aerospace controls engineer a powerful, new 
means of controlling air vehicles. This work explores a new method of designing and 
implementing controllers and guidance systems for autonomous control of air vehicles 
utilizing a GPS integrated guidance, navigation and control system. This is a subject 
of considerable interest when realizing controllers to trank reference trajectories given 
in an inertial reference frame. The design, implementation, and dynamic simulation 
of a precise tracking trajectory controller for an Unmanned Air Vehicle (UAV) is pre- 
sented. This design provides a natural conversion of commands and other measured 
outputs (such as GPS signals) from an inertial reference frame to a body-fixed refer- 
ence frame. This achieves automatic recruiting of the actuators while preserving the 
properties of the original design (linearization principle). 



m 



c,/ 

TABLE OF CONTENTS 

I. INTRODUCTION 1 

II. DEVELOPMENT OF THE DYNAMIC MODEL 4 

A. REFERENCE FRAMES 4 

1. Local Tangent Plane Reference Frame 4 

2. Body-Fixed Reference Frame 5 

3. Flight Path or Wind Reference Frame 6 

B. COORDINATE TRANSFORMATIONS 7 

C. NOTATION 9 

D. RIGID BODY EQUATIONS OF MOTION 10 

1. Linear Motion 10 

2. Angular Rotation 11 

3. External Forces and Moments 12 

4. The State Space Representation 14 

5. Trim and Linearization 19 

III. THE LINEAR QUADRATIC REGULATOR DESIGN 23 

A. LQR OVERVIEW 23 

B. DESIGN REQUIREMENTS 27 

C. THE SYNTHESIS MODEL 28 

D. THE DESIGN PROCESS 32 

E. LQR CONTROLLER PERFORMANCE 34 

IV. CONTROLLER IMPLEMENTATION ON THE NONLINEAR PLANT 41 

A. D-IMPLEMENTATION 41 

B. ^-IMPLEMENTATION OF THE CONTROLLER IN SIMULINK 47 



IV 



MONTEREY CA 93943-510! 



C. GENERATION OF THE TRAJECTORY COMMANDS 48 

D. STATE FEEDBACK TO OUTPUT FEEDBACK 54 

1. Sensor Modeling 54 

2. Kalman Filtering 58 

E. INTEGRATION OF THE FULL NONLINEAR SIMULATION . . 60 

V. APPLICATION TO THE CONTROL OF BLUEBIRD 63 

A. P-IMPLEMENTED CONTROLLER PERFORMANCE CHARAC- 
TERISTICS 63 

B. AN AIRPORT DEPARTURE AND ARRIVAL FLIGHT SIMULA- 
TION 72 

VI. CONCLUSIONS AND RECOMMENDATIONS 78 

A. CONCLUSIONS 78 

B. RECOMMENDATIONS 79 

APPENDIX A: MATLAB FILES 81 

APPENDIX B: SIMULINK FILES 98 

APPENDIX C: P-IMPLEMENTATION PROOFS REFERENCED ... 100 

REFERENCES 106 

INITIAL DISTRIBUTION LIST 107 



LIST OF TABLES 



3.1 EIGENVALUES OF BLUEBIRD 29 

3.2 EIGENVALUES OF BLUEBIRD WITH YAW DAMPER 30 

3.3 TRANSMISSION ZEROS OF SYNTHESIS MODEL 34 

3.4 EIGENVALUES OF THE FEEDBACK SYSTEM 35 

3.5 COMMAND AND CONTROL BANDWIDTHS 38 

4.1 ACCELEROMETER CHARACTERISTICS 55 

4.2 ANGULAR RATE SENSOR CHARACTERISTICS 56 

4.3 INCLINOMETER AND MAGNETOMETER CHARACTERISTICS 56 



vi 



LIST OF FIGURES 



2.1 Local Tangent Plane Coordinate System 5 

2.2 Body-Fixed Coordinate System 6 

2.3 Wind or Flight Path Reference Frame 7 

2.4 SIMULINK Nonlinear Sixteen State Dynamic Model of an Air Vehicle 20 

2.5 SIMULINK Nonlinear Nine State Dynamic Model of an Air Vehicle . 21 

3.1 Standard LQR feedback configuration 24 

3.2 Feedback Configuration for Root Locus Analysis 27 

3.3 Synthesis and Analysis Model 31 

3.4 Control-Loop Bandwidth: Elevator Channel 35 

3.5 Control-Loop Bandwidth: Throttle Channel 36 

3.6 Control-Loop Bandwidth: Aileron Channel 36 

3.7 Comm and- Loop Bandwidth: X Position Channel 37 

3.8 Command- Loop Bandwidth: Y Position Channel 37 

3.9 Command-Loop Bandwidth: Z Position Channel 38 

3.10 Bank Angle and Heading Response to Ramp in Y Command 39 

3.11 Pitch Angle and Altitude Response to Ramp in Z Command 39 

3.12 Trajectory Error Due to a Constant Wind Disturbance 40 

4.1 Block diagram of the nonlinear controller K, (A) 46 

4.2 V- Implementation of controller on Bluebird 49 

4.3 Commanded Trajectory Logic Block 50 

4.4 Example of Commanded Trajectory Revision 52 

4.5 Generation of Trajectory Commands 53 

4.6 Angular Rate Sensor 55 

vii 



4.7 Sensor Models in SIMULINK 57 

4.8 Frequency Response of Position Filter 59 

4.9 Frequency Response of Euler Angle Filter 60 

4.10 SIMULINK Diagram of the Full Nonlinear Simulation 62 

5.1 Test Trajectory #1 65 

5.2 Trajectory #1: Position Error and Euler Angles 66 

5.3 Trajectory #1: Velocity Data 67 

5.4 Trajectory #1: Control Activity 68 

5.5 Test Trajectory #2 69 

5.6 Trajectory #2: Position Error and Euler Angles 70 

5.7 Trajectory #2: Velocity Data 71 

5.8 Trajectory #2: Control Activity 72 

5.9 Trajectory #2: Position Error for Varying Turn Rates 73 

5.10 Trajectory #2: Position Error for Varying Wind Velocities 74 

5.11 Departure and Arrival at an Airfield 75 

5.12 Airfield Scenario: Position Error and Euler Angles 76 

5.13 Airfield Scenario: Control Activity 77 



Vlll 



ACKNOWLEDGMENT 



There axe a few special people I would like mention without whom this project 
would not have been completed. It was a true pleasure to have studied aeronautics 
under Professor Isaac Kaminer and Professor Rick Howard. Their advice and profes- 
sional counsel were invaluable. I am especially grateful to Professor Isaac Kaminer 
whose vision and expertise axe directly responsible for the creation of the superb 
avionics lab where this work was accomplished. Finally, a special thanks to five in- 
dividuals who sacrificed the most. To Nicole, Katie, Eric, and Patty, you were the 
best while understanding the least. And to Patricia, thank you for your exceptional 
patience; without you I would not have been able to complete this project. 



IX 



I. INTRODUCTION 



The advent of GPS has afforded the aerospace controls engineer a powerful new 
means of controlling air vehicles. Current guidance schemes rely in some part on 
ground based radars, navigational aides, beacons, localizer beams, etc. Guidance 
and control of the air vehicle have been developed seperately and then combined in 
a somewhat adhoc process. Precise tracking of inertial fixed trajectories using an 
onboard GPS integrated GNC suite affords a quantum leap in autonomous flight 
capabilities. The greatest impact of the proposed technology is expected in the area 
of trajectory tracking control for autonomous unmanned vehicles, and automatic 
approach and landing of manned vehicles dicussed next. 

Control of the commercial air traffic throughout the country continues to become 
more and more demanding as an increased number of vehicles vie for limited access 
to the major commercial aviation hubs. Sophisticated and expensive ground based 
radar control facilities employ a large number of personnel to individually instruct 
the pilots of these aircraft on the trajectories they are to fly. Precise control of the 
aircraft trajectory such as required by an approach into an airfield requires constant 
attention from a ground based air traffic controller. Furthermore, ATC ability to 
effectively control the aircraft is influenced by ground based radar coverage and navaid 
equipment available and is often negatively influenced by atmospheric conditions. 
Airfields with limited resources are often unable to take-off and land aircraft requiring 
instrument departures and arrivals. 

Flight patterns around major aviation hubs are, in general, inertial based tra- 
jectories. That is, an aircraft is required to track a certain path over the ground while 
adhering to a certian altitude schedule, irrespective of air mass disturbances. In some 
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cases, such as on final approach to land or when the aerodrome is situated among 
significant terrain or cultural development, precise adherence to the desired trajectory 
is crucial for flight safety. In any case, commanding an inertial trajectory directly 
utilizing a GPS integrated GNC system could prove to be more cost effective and 
accurate than current methods. Furthermore, such a guidance scheme could open up 
many more airports to significant commercial air traffic without requiring the capital 
investment and maintenance of ground based radar facilities. 

Unmanned air vehicles can be a cost effective means of power projection. Addi- 
tionally, in some cases human physiological limits may prove to be the limiting factor 
in the performance of an air vehicle. The precision delivery of munitions becomes 
of paramount importance as weapons and weapon delivery platforms continue to in- 
crease in cost, thus limiting their numbers. All of these concerns can be addressed 
by autonomous air vehicles utilizing a GPS aided guidance, navigation and control 
suite. 

All of these applications have a common thread r unning through them. While 
the particulars of the vehicles may vary considerably, the intent is to acheive au- 
tonomous control of their trajectory. As a proof of concept, this work presents a new 
design process for the synthesis of a guidance, navigation, and control system for a 
UAV named Bluebird. The function of the this GNC system is to track inertial tra- 
jectories. Bluebird is a UAV operated at the Unmanned Air Vehicle Lab at the Naval 
Postgraduate School. It has a 12.5 foot wingspan and a 20 pound payload capability, 
and is currently being equipped with a full avionics suite, including IMU, GPS, and 
air data sensors. 

The design process began with the development of a nonlinear dynamic model of 
Bluebird implemented in SIMULINK. A typical cruise flight condition was chosen as 
the point for linearization. After linearization of the nonlinear model, the work cen- 
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tered on the design of a linear controller. LQR (Linear Quadratic Regulator) synthesis 
approach was used since it provides an intuitive means of synthesizing a multivariable 
controller within the framework of real world design constraints. Following the de- 
sign of the LQR controller, the challenge of implementing the linear controller on the 
nonlinear plant was addressed. A novel method of converting commands and outputs 
from inertial to body reference cooridinates was used [Ref. 10]. This method achieves 
automatic recruiting of the actuators, while preserving a certain linearization prop- 
erty. Next, the accuracy of the nonlinear simulation was enhanced with the addition 
of high fidelity models of sensors used onboard Bluebird. Additionally, Kalman filters 
were designed in order to provide optimal state estimates. Finally, the performance 
of the controller was evaluated in simulations with the full nonlinear model. 
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II. DEVELOPMENT OF THE DYNAMIC 

MODEL 



The development of an integrated guidance, navigation, and control system 
required a high fidelity nonlinear model of the aircraft dynamics. The discussion 
begins with an explanation of nomenclature, abbreviations, and a definition of frames 
of reference. 

A. REFERENCE FRAMES 

Three different reference frames are used in this report. They are: 

• Local Tangent Plane or Inertial Reference Frame 

• Body-Fixed Reference Frame 

• Wind or Flight Path Reference Frame 

1. Local Tangent Plane Reference Frame 

The position of the air vehicle must be maintained with respect to the local 
tangent plane coordinate system. This coordinate system is formed by extending a 
ray from the center of the earth to its surface. A plane is attached tangent to the point 
of intersection of the ray with the Earth’s surface. While it is somewhat arbitrary, for 
our purposes here it will be convenient to define the positive x direction as pointing 
east, the positive y direction as pointing north, and the positive z direction as pointing 
up. This is depicted in Figure 2.1. 

For the purposes of this development, the rotation of the earth and its 
associated Coriolis’ forces can be ignored and the local tangent plane reference frame 
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can be considered to be an inertial reference frame. In this work, {/} is used to 
represent the inertial reference frame. 

2. Body-Fixed Reference Frame 

The body-fixed reference frame is a right hand orthogonal system with 
the origin at the center of gravity of the air vehicle. The positive x direction points 
towards the nose. The positive y direction points out the right wing and the positive z 
direction points towards the bottom of the air vehicle. The velocity of the air vehicle 
with respect to the inertial reference frame, resolved along the x, y, and z axis of the 
body-fixed reference frame, are termed u, v, and w, respectively. The angular rate of 
rotation of the air vehicle with respect to the inertial reference frame, resolved in the 
body-fixed reference frame, are called p, q, and r, respectively. Positive values for the 
forces, moments, angular rates, and linear velocities in the body-fixed reference frame 
are shown in Figure 2.2. The abrieviation, {B}, is used to represent the body-fixed 
reference frame. 
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Figure 2.2: Body-Fixed Coordinate System [Ref. 11] 

3. Flight Path or Wind Reference Frame 

The wind reference frame is also a right hand orthogonal system with its 
origin at the center of gravity, c.g., of the air vehicle. The x axis is aligned with the 
velocity vector of the air vehicle. The orientation of the wind reference frame with 
respect to the body-fixed reference frame is defined in terms of the angles a and (3. 
The equations for a and (3 are given below. 



and 



a = tan i (w/u ) 


(2.1) 


(3 = sin~ l (v/V) 


(2.2) 
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where the vectors u, u, w, and V are velocity components of the air vehicle defined in 
Figure 2.3. The abrieviation, {W}, is used to represent the wind reference frame. 




Figure 2.3: Wind or Flight Path Reference Frame [Ref. 11] 

B. COORDINATE TRANSFORMATIONS 

In order to use these three coordinate systems, one must be able to transform 
between them freely. The Euler angles, 0, and 'F, termed roll, pitch, and yaw, 
are defined in order to express the orientation of the body-fixed reference frame with 
respect to the inertial reference frame. For the purposes of this development, a 3-2-1 
Euler angle transformation will suffice as its singularity occurs at 0 equal to 90 de- 
grees. The 3-2-1 transformation is given without explanation but a good development 
of Euler angle transformations in general can be found in [Ref. 6]. The nature of 
the angular rotation is more apparent when the transformation is expressed as the 
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product of three rotation matricies. In the case of a 3-2-1 rotation sequence, the three 
matrices in Equation 2.3 correspond to rotations about the yaw, pitch, and roll axes 
of the air vehicle. Of course, the three matrices can be multiplied out for an analytic 
result contained in a single matrix , although the resulting matrix is somewhat busy 
to inspect. In any case, the transformation between a free vector resolved in the in- 
ertial reference frame and the same vector resolved in the body-fixed reference frame 
is given by: 



r V = 



cos '5 


sin 


0 ‘ 




cos 0 


0 


— sin0 




' 1 


0 


0 ‘ 


— sin $ 


cos $ 


0 




0 


1 


0 




0 


cos $ 


sin $ 


0 


0 


1 




sin0 


0 


cos 0 




0 


— sin$ 


cos $ 



B 



V (2.3) 



where l V is a free vector resolved in {/} and B V is the same vector resolved in {5}. 
The inverse is also defined. Conveniently, since the transformation is orthonormal, 
the inverse is simply the transpose of the rotation matrices shown in Equation 2.3. 

Not all, transformations, however, axe orthonormal. Of particular interest is 
the case of angular rotation rates. The body-fixed reference frame’s angular rate of 
rotation with respect to the inertial reference frame can be related to the rate of 
change of the Euler angles by a transformation matrix. The development is straight 
forward and is fully explained in [Ref. 11]. The final transformation matrix from p, 
q, r to the time rate of change of the Euler angles, $, 0, is given by: 



' $ ' 




' 1 


sin $ tan 0 


cos $ tan 0 


0 


= 


0 


cos $ 


— sin$ 


. * . 




0 


sin$ sec0 


cos $ sec 0 



P 

? • 
r 



(2.4) 



By integrating Equation 2.4, the time history of the Euler angles can be obtained. 



Aerodynamic forces and moments axe often calculated using stability and control 
derivatives defined with respect to the wind reference frame. The angles, a and /3, 
define the orientation of the wind reference frame to the body-fixed reference frame. 
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Therefore, a transformation matrix can be obtained that relates a free vector, such as 
lift or drag, resolved in {V^} to the same vector resolved in {£}. The transformation 
is expressed as: 



B V = 



cos a cos /? — cos a sin /? — sin a 
sin j3 cos /3 0 

sin a cos /3 — sin a sin /3 cos a 



Wy 



(2.5) 



where B V is a free vector resolved in { B } and W V 



is the same vector resolved in 



{W}. 



C. NOTATION 

Some standardized abbreviations will simplify the development of the nonlinear 
kinematic model of the air vehicle. This short-hand is used in the field of robotics 
where multiple frames of reference are co mm on [Ref. 5]. 

• I P cg represents the position vector from the origin of the local tangent plane to 
the center of gravity of the air vehicle. 



• B Vcg and B acg represent the velocity and acceleration, measured at the center of 
gravity of the air vehicle, with respect to {/}, resolved in {5}. The components 
of B v cg are commonly termed u, v, and w. 



• ^ and l a cg represent the velocity and acceleration, measured at the center of 
gravity of the air vehicle, with respect to {/}, resolved in {/}. 

• b u>b is the angular velocity of the {f?} coordinate system with respect to {/}, 
resolved in {B}. The components of b u>b are commonly termed p, q, and r. 

• represents the angular velocity of the {B} coordinate system with respect 
to {/}, resolved in {/}. 
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• gR represents the transformation matrix used to express a free vector resolved 
in {B}, in {/}. The inverse is represented by f R 

• fyR represents the transformation matrix used to express a free vector resolved 
in {W}, in {5}. The inverse is represented by 'q R. 

• B F and B N denote the total external inertial force and moment acting on the 
body resolved in {B}. 

• 1 F and 1 N denote the total external inertial force and moment acting on the 
body resolved in {/}. 

• B L is the inertial angular momentum of the body resolved in {5}. 

• 1 L is the inertial angular momentum of the body resolved in {/}. 

• Given a vector v, its derivative with respect to {B} is denoted as £(v) 
and its derivative with respect to {/} is denoted as (u) 

D. RIGID BODY EQUATIONS OF MOTION 

In general, an avionics suite on a modem air vehicle utilizes a strapdown IMU. 
A strapdown IMU, as the name implies, maintains a constant orientation in the body- 
fixed reference frame. The output of the sensors on the IMU are resolved in {5}. 
Therefore, among other reasons, it is most convenient to develop the equations of 
motion of the air vehicle in the body-fixed reference frame. 

1. Linear Motion 

An application of Newton’s Law to linear motion of a body states that 
the total external force applied to a body is equal to the mass of the body times its 
inertial acceleration. This could be written in the inertial reference frame as: 
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where 



r F 



m r a, 



i 



a = 




or in the body-fixed reference frame as follows: 



( 2 . 6 ) 



b F = m B a 

= m f&Vcg 

= m B v cg . (2.7) 

Coriolis’ theorem can be used to relate the inertial and body accelerations of the air 
vehicle as follows: 



B V cg = -fo Bv cg + B UB X B V eg , (2.8) 

where the difference in the derivatives is explained in Chapter II, part C. Equation 2.8 
can be substitued into Equation 2.7 in order to obtain the desired expression for the 
sum of the external inertial forces resolved in the body-fixed reference frame. 



B F = m (■ j B Vcg + B ub x B v cg ) 

= m ^ B Vcg + m ( b u b x B v cg ). 



(2.9) 



2. Angular Rotation 

Euler’s law for the conservation of angular momentum at the center of 
gravity states that: 
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( 2 . 10 ) 



l L — l N 

lJcg “ iy cgi 

where 1 L cg is the angular momentum of the air vehicle with respect to {/} and 1 N cg 
is the total moment applied to the air vehicle. Equation 2.10 can be written in the 
body-fixed reference frame as: 



B L<, = ffl'iV 

Coriolis’ theorem can be used again to expand B L cg obtaining: 



( 2 . 11 ) 



B L cg = ± B Lc, + b u b x B L cg . (2.12) 

It can be shown that the angular momentum, B L cg , of the air vehicle is the product 
of an inertia tensor, defined as J B , and the body’s angular velocity, B u> B , where we 
ignored all spining elements. Substituting this definition of B L cg into Equation 2.12, 
results in: 



B L cg = ~^(J b b u b ) + B IjO B x J b b u b . (2.13) 

Recall that B = f R 1 N cg = B N^. Using this relationship, Equation 2.13 can be 
equivalently expressed as: 



B N cg = -j-(J b b w b ) + b u b x J b b u b (2.14) 

at 

3. External Forces and Moments 

Equation 2.9 and Equation 2.14 from the preceding sections can be com- 
pactly expressed as follows. 
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By bringing derivative terms in Equation 2.15 to the left hand side, we obtain, 



d_ 

dt 



B V C g 




- b vb x 


B U) B 




“ Jq 1 X Jb B Ub + j B B N 



(2.16) 



Next, the forces and moments in Equation 2.16 can be expanded as follows: 



■ Bp ' 




B FgraVITATIONAL + B FpROPULSIVE + B F AERODYNAMIC 


b n 




B Npropulsive + b N aerodynamic 



(2.17) 



where 



B F, GRAVITATIONAL 
B FpRopuLsrvE 
B F AERODYNAMIC 
B N PROPULSIVE 
B N AERODYNAMIC 



force due to gravity 
force due to engine's thrust 
force generated by aerodynamic surfaces 
moment generated by engine's thrust 
moment generated by aerodynamic surfaces 



The gravitational force expressed in {/} is given by: 

! F gravity = 

where ”g” is the gravitational constant. Then 



0 

0 

mg 



B Fgravity = f R 1 Fgravity- (2-18) 



The expansion of the propulsive forces and moments is simplified by consid- 
ering the case of centerline thrust. In that case, no external moments are generated, 
i.e., B Npropulsive = 0, and propulsive forces can be expressed as: 



' T ‘ 

b Fprop = 0 , 

0 



(2.19) 



where T represents the thrust of the powerplant. 

Aerodynamic forces and moments are commonly calculated using nondi- 
mensional stability and control derivatives. These derivatives are obtained by approx- 
imating the aerodynamic forces and moments acting on the air vehicle using a Taylor 
series expansion about a given trim point. Typically, values for these derivatives are 
available for the first order terms of the expansion only. Sometimes, a few second 
order terms are available, such as the terms associated with a and /3 [Ref. 7]. All 
other higher order terms are usually ignored in this approximation. 

In general, the aerodynamic forces and moments acting on the air vehicle 
are computed as follows. The nondimensional stability and control derivatives are 
dimensionalized by multiplication by the appropriate constants, such as wing span, 
dynamic pressure, chord, etc. The dimensional derivative is then multiplied by the 
perturbations of each aerodynamic variable or control deflection from its nominal 
trim point. The summation of the forces and moments due to all of the aerodynamic 
variables and control deflections, in addition to the trim value of the forces and 
moments, results in the total aerodynamic force and moment acting on the air vehicle. 

4. The State Space Representation 

In order to implement the dynamic model in a state-space form suitable for 
numerical simulation, the states of the model need to be chosen. This is somewhat 
arbitrary and many choices will work so long as consistency is maintained in the 
approach. 

As was evident from the development of the rigid body equations of motion, 
the body-fixed reference frame is the most convenient coordinate system in which to 
define the states. The first three states are defined as the inertial velocity of the air 
vehicle resolved in the body-fixed reference frame. These are abrieviated as u, v, and 
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w or more compactly as B v cg . The fourth through sixth states axe defined as the 
angular velocity of the air vehicle with respect to {/} resolved in {B}. These are 
abrieviated as p, q, and r or more compactly as b <jJb- 
Control inputs axe represented by the vector A: 



A = [6 e ,S r ,6 a ] (2.20) 

where S e , S r , and S a axe the elevator, rudder, and aileron inputs, respectively. 
Control inputs for the throttle axe represented by S tr t- 

Typically, the terms in the Taylor series expansion for the aerodynamic 
stability derivatives axe partial derivatives with respect to u/U, a, /?, p, q, and r, 
where U is the magnitude of the air vehicle’s velocity vector and a and /? axe the 
angles defining the orientation of {B} with respect to {W} [Ref. 3]. The last three 
variables, p, q , and r, axe states of the model. The first three can be represented as 
a combination of states in the model as follows. First note that, 



* U ’ 
0 
0 



b K v, 



cgi 



( 2 . 21 ) 



and for small values of angles a and /?, a is approximately equal to w/U and /? is 



approximately equal to v/U. 

It turns out, the stability derivatives can be placed in matrix form as 

follows: 



Clu 


C L0 




Cl, 


Cl, 


c Lr 


Cyu 


C Yfi 


Cy„ 


Cy, 


Cy, 


Cy t 


Cdu 


Cd 0 


Cd* 


Cd, 




c Dr 


Clu 


Oh 


c, a 


Ci, 


Ci, 


C lr 


Cm u 


Cr 


c ma 


Cm p 


Crriq 


C m r 


cw 


Cn, 


c na 


Cn? 


C nq 


C„ r 



Similaxily, for the control derivatives: 
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dC_ 

dA 



J 





Cl (t 


C L, a 


C Y" 


Cy (t 


C Y" 


C D St 


Cd (t 


C D, a 


Cl" 


Ci Sr 


Ci (a 


Cm St 


Cm fr 


Cm fa 


Cn fe 


Cn Sr 


Cn ta J 



where x' is the vector composed of u/U , a, f3, p, q, and r. Now the aerodynamic 
forces and moments can be expressed as follows: 



B Faero 
B N aero 



= qS 









f „ dC w/ dC lV/ . ac A1 
+ a? M 1 + x + 



dx' 



dA 



( 2 . 22 ) 



where A/', q, and S axe matrices used to dimensionalize the stability and control 
derivatives and convert the state vector, x, to x': 



x 




u 



v w p q 




q = dynamic pressure, 



S = diag{— 5, S, -S, Sb, Sc, Sb}, 



x' = M'x, 



M' = diagfl/Fr, 1/V T , l/V T , b/2V T , c/2V T , b/2V T }, 
x' — M'x, 

M' = diag{0, c/(2Vt), b/ (2Vr)»0, 0, 0}. 



Equation 2.16 can now be further expanded using expressions of the forces 
and moments derived above. 
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d_ 

dt 



where 



b , 

B 



'eg 

UB 



- B U B X 


0 




0 


- b Jb\ b u b x 1 


’ Jb b 

■ i_ 
m 






0 


s F 1 


j \ b Fgrav 


+ 


b n\ 


U 0 . 



B J~b 1 



u cg 

Ub 
B p 

B N 



+ 



$trt+ 



(2.23) 



• qS { C F0 + fgM'x + §M'i + IfA } | | . (2.24) 

Notice that there is a state derivative term on the right hand side of Equa- 
tion 2.23 due to the second order terms in the Taylor series expansion of the aerody- 
namic forces and moments in Equation 2.24. By bringing it to the left hand side of 
Equation 2.24 and combining terms, we obtain: 



0 

0 &R 



£ 

dt 



B, 



J cg 

3 U B 



-rill 



- B u B x 



0 - B J B \ B »BX B JfuB) 



Mr'&TiS&M' ] [ ll 



eg 

Ub 



+ Mj 



■U 



b Fgrav 

0 



b Fprop 

0 



6 trt + &TqS(C Fo + SfeA) 



}}■ 



(2.25) 



where: 



B rp 

w 1 — 



&R 0 
0 &R 



and Mi = 



m 0 

0 b Jb 



and 



X — h~ MT'ZrTqS^M'. 



(2.26) 



Equation 2.25 expresses the derivative of the first six states of the nonlinear model 
in matrix form. It is solved in the user defined MATLAB Fen block, state-deriv.m , 
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which is available for inspection in Appendix A. The physical parameters specific 
to Bluebird are stored in a MATLAB .m file, Bluebird-data. m, and are called from 
within the function, state-deriv.m . In order to change the vehicle being modeled, 
one simply changes the constants in Bluebird-data.m . 

Next, the Euler angles were added as the additional three states of the 
model. The time history of the Euler angles is defined in Equation 2.4 and written 
in compact form as: 



A = S( A)Vb, 



(2.27) 



where 



and 



S{ A) = 



A = 



$ 

0 



1 sin $ tan 0 cos $ tan 0 
0 cos $ — sin $ 

0 sin $ sec 0 cos $ sec 0 



The user defined MATLAB Fen block, eu/.m, solves Equation 2.27, the details of 
which can be found in Appendix A. 

Finally, the inertial position of Bluebird can be computed as follows. 



'Pc =' (2.28) 

The rotation matrix, B R, is implemented in a MATLAB Fen block, 
posbSi.m, in order to convert the vector B v cg to the vector I v cg . The vector 
is then integrated to obtain a time history of the position of the air vehicle in the 
local tangent plane. To increase fidelity of the simulation, a first order model of the 
four actuators, b elevator , ^rutWcr) ^aiicroni b tht ottic i is included with a time constant 
of 1/12. The resulting nonlinear dynamic model now has sixteen states which are 
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computed using Equations 2.25, 2.27, and 2.28; summarized here for clarity: 
states = [izuu;p 9 r$0<P Xpos Ypos Zpos S Cpot 8 Tpo , 6 apo . 6 tpot 

= [X bjt A t p t A<.] T , 



£ 

dt 



v c g ] _ -1 / [ - B U>B* 0 

u) B * 1 o 

Mp&TqS&M' j 



B 



- b J b 1 ( b u b x b J b b uj b ) 

b Fgrav 

0 



+ 



+ UT l | 



B FpBOP 

0 



Ub 

<« + &r55(C F „ + ^A) 



)}■ 



(2.29) 



A = S(A) b u b , 



(2.30) 



J Ai = (2.31) 

The SIMULINK diagram of the nonlinear model is shown in Figure 2.4. 

5. Trim and Linearization 

For linear controller design, the nonlinear equations above must be trimmed 
and linearized for a typical cruise flight condition for Bluebird. A SIMULINK tool is 
available which can be used to find equilibrium points of nonlinear dynamic models. 
The user specifies which states and control inputs are to remain fixed along with 
their stationary values and the trim routine searches for values of the state and input 
vector for which the derivative of the state vector equals zero. With the trim condition 
known, another SIMULINK tool perturbs the states around the specified trim point 
in order to find the rate of change of the states and control inputs (Jacobian). The 
resulting linear model is returned in state space format. Since Equation 2.28 can only 
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Figure 2.4: SIMULINK Nonlinear Sixteen State Dynamic Model of an Air 
Vehicle 

be trimmed for v = 0, not a typical flight condition, we will trim Equations 2.25 
and 2.27, and then include Equation 2.28 for linearization. 

We are interested in trimming the model in velocity. Hopefully, the deriva- 
tive of the position states will never equal zero in flight. Also, in trim, the control 
inputs are equal to the actuator positions. Therefore, actuator states are also removed 
from the nonlinear model. The nine state nonlinear model of Bluebird used for trim 
is shown in Figure 2.5. 

A typical cruise flight condition for Bluebird is given by: 
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Mux 




1 




s 












Figure 2.5: SIMULINK Nonlinear Nine State Dynamic Model of an Air 
Vehicle 

• flight speed equal to 73 feet per second 

• flight path angle equal to zero 

• wings level attitude 

The nonlinear model depicted in Figure 2.4 was trimmed at this condition. 
The trim values of the nine states, u, v, w, p, q, r, $, 0, and tf, and the four 
control inputs, 6 e , 6 r , S a , and 6 t were returned. While the sixteen state nonlinear 
model cannot be trimmed in position, it can be linearized at an arbitrary position. 
The origin of {/} is conveniently chosen. These values were then used to linearize 
the complete nonlinear model, including position and actuator states. The resulting 
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linear model of Bluebird and numerical values for the trim condition are included in 
the Appendix B. 
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III. THE LINEAR QUADRATIC REGULATOR 

DESIGN 



The previous chapter developed a nonlinear model of Bluebird, which was used 
to derive a linear model for a cruise flight condition. In this chapter, a linear dynamic 
controller is developed to provide trajectory tracking for the linear model. LQR 
methodology was selected to design the controller. Based on design requirements, an 
intuitive means of manipulating the LQR gains is presented. See [Ref. 9] for details. 
The following is a brief review of the properties of an LQR controller utilized in this 
design process. 



A. LQR OVERVIEW 

Consider the linear system 



Q = 



X 

i z 



y 



Ax + Bu 
C\x + D\u 
x 



where x 6 /?", it € BT 1 and z € BP. 

Suppose C\D\ = 0 and D\ is full column rank. Then, 



z T z = x T C*C\x + u t D[D\U. 



Define a cost, J, as follows: 



(3.1) 



J— J(z T z)dt = J (x T C^C'ix + u T D[ D\u)dt (3.2) 

and let Q = C?C U and R = Df Lh. Note: Q >0 and R > 0. 

Assume ( C\,A ) is observable and (A, B) is controllable. Consider Figure 3.1. 
The standard LQR problem is to find a controller, u = K(s)x , such that the feed 
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i 



back system in Figure 3.1 is internally stable and J is minimized. It turns out, one 
such controller uses a constant gain, u = —Kx, where 



K = - R~ 1 B t P , (3.3) 

and P solves the Alegbraic Riccatti Equation: 

A t P + PA- PBR~ 1 B t P + Q = 0. (3.4) 

Figure 3.1 shows the feedback interconnection of the plant Q and the controller K. 
Here the inputs w can include commands and disturbances. 



Q ci 




Figure 3.1: Standard LQR feedback configuration. 



It turns out that the controller, K, has guaranteed simultaneous phase and gain 
margins of no less than 60 degrees and 3dB, respectively [Ref. 12]. Furthermore, the 
controller has asymptotic properties which are exploited in the design process and 
are discussed next. 

Define the Hamiltonian matrix, H, as: 



H = 



A -BR~ l B T 
Q -A T 



Let T be given by: 
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where P solves the Riccatti equation, Equation 3.4. 
Note, 



T-1 UT - r A - BR~ 1 B T P -BR~ 1 B t P 

[ 0 -A t + PBR~ 1 B t 

Therefore the eigenvalues of H are the roots of the following polynomials: 

det(sI-H) = det(sl —T~ X HT) = det(sI-A+BRr 1 B T P)det((sI+A T -PBR~ 1 B T ) 

Clearly, the eigenvalues of the Hamiltonian consist of the eigenvalues of Q and their 
unstable reflections about the imaginary axis. Let R = pR\, R\ > 0, then; 



si -A (l/p)^ 1 ^ 
0 si A A 



det(sl -H) = 

It can be shown that the det(sl — H) can be equivalently expresses as [Ref. 9]: 



det(sI-H) = - rdet(sI-A)det(-sI-A)det(I+(l/p)Ri 1 B T {-sI-A T )- 1 CTC 1 (sI-A)- l B ) 

Let 

^(s) = det(sl — A) 



and 



6(s) = Ci{sl - A)~ l B. 



Then 



detfsl - H) = -r<p(s)^(-s)iet(I + (l/p)Rp6(-s)6(s))). 



(3.5) 



The SISO example will best demonstrate what occurs to the eigenvalues of Q d 
as p is varied from 0 to oo. Let 



0(a) = ip(s)/<j>(s) 

where 0(s) axe the zeros of 0(s). Then, 

det(sl -H) = —l n <f>(s)<f>(—s)(I + (l/p)0(s)0(-s)/^(s)^(-s)). 

It follows that the eigenvalues of H are the roots of, 

<j>(s)<f>(-s) + (1 / p)rJ>{3)ij>{-3). (3.6) 

Consider the feedback system shown in Figure 3.2. This is standard configuration for 
SISO root locus analysis and its characteristic equation is: 

<l>(s)<f>(-s) + (l/p)rJ>{s)rJ>(-3), 

exactly the same as Equation 3.6. Since we know that the stable eigenvalues of H 
are the eigenvalues of Q standard root locus techniques show that 
as p goes to 0 [Ref. 9]: 

• p eigenvalues of Q ^ go to the stable zeros of Ci(sl — A)~ l B or the stable 
reflection of the unstable zeros of Ci(sl — A)~'B, where p is the number of 
zeros of C\(sl — A)~ l B. 

• n-p eigenvalues of Q d go to -00 in Butterworth patterns, 
as p goes to 00 

• n eigenvalues of Q d go to the stable eigenvalues of A and the stable reflections 
of the unstable eigenvalues of A. 
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i/p 






X 








ip(-s)xj>(s) 







Figure 3.2: Feedback Configuration for Root Locus Analysis. 

B. DESIGN REQUIREMENTS 

These properties of an LQR controller were utilized in a design process in order 
to meet the following design requirements. 

1. Zero Steady State Error 



• Achieve zero steady state values for all error variables in response to ramp 
commands in position along the x, y, and z inertial axes. Note that a ramp 
command in postion corresponds to a constant heading, constant velocity 
trajectory. 

2. Bandwidth Requirements 

• The input-output command response bandwidth (command-loop band- 
width) along any of the three command channels should be no greater 
than 1 radian per second and no less than 1/10 radian per second. 

• The control-loop bandwidth should not exceed 12 radians per second for 
the elevator and aileron actuators, and 5 radians per second for the throttle 
actuator. These numbers represent 80 % of the corresponding actuator 
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bandwidths and shall ensure the actuators are not driven beyond their 
linear operating range. 

3. Closed Loop Damping 

• The dominant closed loop eigenvalues should have a damping ratio of at 
least 0.7. 

C. THE SYNTHESIS MODEL 

The synthesis model is the primary interface between the control design 
and the LQR algorithm. At the heart of the synthesis model is a linear model 
of Bluebird developed in Chapter II. 

Bluebird has four control inputs, namely elevator, rudder, ailerons, and 
throttle. The elevator and throttle are natural choices for controlling x and 
z position in steady state. The remaining two control inputs could be used 
to control the lateral variable (y position). Both rudder and aileron provide 
means of generating accelerations in the lateral plane. In fact, rudder is more 
effective at generating sideslip than aileron. In the linear plant, lateral position 
is the double integral of lateral acceleration. Subsequently, the resulting LQR 
controller will attempt to use rudder to null out errors in lateral position, i.e., 
to turn the plane. However, the desired controller response is to bank to turn 
using ailerons and to use rudder for turn coordination. Furthermore, in the 
presence of wind, it is desired that Bluebird fly wings level, crabbed into the 
wind, rather than use a wing down, top rudder technique. For these reasons, 
the rudder was removed as a control input to the linear model. 

As can be seen from Table 3.1, the dutch roll mode of Bluebird is lightly 
damped. This light damping of 0.111 could pose a performance problem. Since 
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rudder is available and not used in the design of the trajectory controller, the 
nonlinear model of Bluebird was modified to include a yaw damper for improved 
dutch roll damping. Yaw rate was fed back to the rudder through a constant 
gain block with a value of 0.55. Additionally, the rudder was removed as an 
external input. Note that Bluebird is still fully controllable with the remaining 
three control inputs. 

TABLE 3.1: EIGENVALUES OF BLUEBIRD 



Mode 


Frequency 


Damping 


Longitudinal 


rad/sec 




Short Period 


5.9 


0.735 


Phugoid 


0.497 


0.0344 


Lateral-Directional 






Dutch Roll 


2.4 


0.111 


Spiral 


0.0384 


-1 


Roll Response 


-4.572 


1 



The no nlin ear model of Bluebird with three inputs and integral yaw damper 
was linearized, as per Chapter II, returning the linear model, 



- f x = Ax + Bu 

e =\v = X 

where 

x = ^ u v w p q r $ Q W Xpos Ypos Zpos S epo , 6 apo , 6 tpo , ] 
and 

U — | 6 elevator & aileron & throttle J. 



(3.7) 



This linear model was used in the LQR design. The eigenvalues of Bluebird 
with a yaw damper axe given in Table 3.2 where it can be seen that the dutch 
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roil mode has been damped out. Note that the number of states decreases to 
fifteen since the rudder actuator state was removed. 



TABLE 3.2: EIGENVALUES OF BLUEBIRD WITH YAW DAMPER 



Mode 


Frequency 


Damping 


Longitudinal 


rad/sec 




Short Period 


5.9 


0.735 


Phugoid 


0.497 


0.0344 


Lateral-Directional 






Dutch Roll 


2.35 


0.5 


Spiral 


0.1788 


1 


Roll Response 


-4.5686 


1 



Consider Figure 3.3. Here K is the controller to be designed, G is the 
linear model of Bluebird, and the block, 5, within the dotted line is the synthesis 
model. 

The signal, in, represents the commanded trajectory inputs: 



^ J XpOScmd YpOScmJ ZpOScmd ^ cmd ® cmd j 



The signal x\ represents the linear and angular position states in the linear 
model. 



X! 



Xpos Ypos Zpos $ 0 'F 



r 



The signal e represents the errors between the commanded and current trajec- 
tory. The signal z is comprised of the outputs of the matrices, Q , and R. Since 
zero steady state error is desired while tracking a ramp command in inertial 
position, two integrators were placed on each error signal. This also ensures 
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s 




Figure 3.3: Synthesis and Analysis Model 

perfect tracking of constant heading trajectories in the presence of a constant 
wind disturbance. Thus Q was chosen to be; 





‘ 91 


0 


0 ' 




£11 

1 


£12. 

» 




Q = 


0 


92 


0 




£21 

1 


£ 22 . 

9 






. 0 


0 


93 . 




£21 
L l 


£ 22 . 

9 


¥ . 



The values of c** were chosen to place six transmission zeros from u to z at 
appropriate locations. If they axe well chosen, six poles in the closed loop plant 
will move very near to the placed transmisssion zeros. With that in mind, the 
transmission zeros were chosen as appropriate target locations for the poles 
added by the addition of the error states. 

The q xx weightings are used as a mechanism for obtaining the desired 
command bandwidth. Increasing the value of q xx increases the relative propor- 
tion of that error state in the regulated output vector z. The resultant LQR 
gain increases the command bandwidth in that channel in order to move the 
controlled state to its commanded value more quickly. 
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The matrix R is a constant diagonal matrix required to be full rank. Since 
the plant Q has three control inputs, R is of the form, 

( r n 0 0\ 

R = I 0 T22 0 

\ 0 0 r 33 ) 

The elements of R are used as a mechanism for selecting the control bandwidth. 
An increase in r xx increases the relative proportion of that actuator energy in the 
regulated output z. The resultant LQR gain decreases the control bandwidth 
of that control input. 

D. THE DESIGN PROCESS 

Design requirements given in the previous section are SISO in nature. 
They are expressed as bandwidth limitations of the individual actuators and 
rise time and damping characteristics along the command channels. Note, the 
rise time is inversely proportional to the command bandwidth. The following 
LQR design process provided a means of obtaining a multivariable solution to 
achieve SISO design specifications. 

With an appropriate linear representation of Bluebird and a synthesis 
model that incorporated well placed transmission zeros, the design ” knobs” 
were adjusted in order to meet performance requirements. The design "knobs” 
axe the elemental weightings, q xx and r xx , in the Q and R matrices. The design 
process iterated through the following steps. 

1. Initially let q xx equal 1. Iteratively determine weights for R to sat- 
isfy control loop bandwidth requirements. Increasing r xx decreases the control 
bandwidth along that channel. 

2. With R from step 1, iteratively determine weights for q xx to satisfy 
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command loop bandwidth requirements. Increasing q xx increases the command 
bandwidth along that channel and decreases the rise time. 

3. If it is required to increase the damping of a lightly damped mode, use 
an eigenvector decomposition to determine the primary states affecting that 
mode. Include a weighting on the derivative of those states in the output z. 

4. Ensure that control-loop bandwidths axe still satisfactory with the 
values of q xx in step 2. It is possible that all of the performance requirements 
axe not acheivable within control bandwidth limitations. 



5. Connect the LQR controller to the linear plant and evaluate the per- 
formance in terms of command response and disturbance rejection. 

6. Confirm satisfaction of other design requirements, including damping. 

7. If any step is unsatisfactory, go back to the synthesis model and make 
appropriate changes. Transmission zeros may need to be added, moved, or 
deleted. Synthesis model outputs may need to be reevaluated. 



Bode plots were used to determine compliance with the requirements. 
After five iterations through the seven step process, the following values for Q 
and R matrices resulted in a controller design that met design requirements. 





' 0.5 


0 


0 ' 




Q = 


0 


1 


0 






0 


0 


1 





o.oi 

0.(562 



;25 



1 02 
8 

l 2d 
1 0% 0.(5625 

1 » 



R = 



5000 0 0 

0 1000 0 
0 0 1000 



The transmission zeros created in the sythesis model are shown in Ta- 
ble 3.3. 
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TABLE 3.3: TRANSMISSION ZEROS OF SYNTHESIS MODEL 



Channel 


Cxi 


Cx 2 


Cx3 


Freqency (rad/ sec) 


Damping 


Xpo, 


1 


0.2 


0.01 


0.1 


1 


Y p o. 


1 


0.4 


0.0625 


0.25 


0.8 


Zp Oi 


1 


0.4 


0.0625 


0.25 


0.8 



E. LQR CONTROLLER PERFORMANCE 

The eigenvalues of the feedback interconnection of the plant and controller 
axe given in Table 3.4. It is apparent that the zeros created in the synthesis 
model were well placed and attracted the integrators created by the addition 
of the error states. There axe two sets of lightly damped poles. These do not 
present a problem because their frequency is an order of magnitude greater 
than the frequency of the eigenvalues associated with the trajectory commands. 
Notice that the actuator poles did not change, indicating that the control band- 
widths were slower than the actuator bandwidths. Actuator models provide a 
simple means of determining if the control bandwidths exceeded the actuator 
bandwidths. 

Figures 3.4 through 3.6 depict the control-loop bandwidths for the eleva- 
tor, aileron, and throttle. The cross coupling between longitudinal and lateral 
flight controls was so slight that it is not shown due to scale. 

Figures 3.7 through 3.9 depict the command- loop bandwidth for step com- 
mands in inertial position. Notice that there is some coupling between X com- 
mand and Z response; the rest axe essentially uncoupled. 

A summary of the resulting command and control bandwidths achieved 
is presented in Table 3.5. 
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TABLE 3.4: EIGENVALUES OF THE FEEDBACK SYSTEM 



Mode 


Frequency (rad/ sec) 


Damping 


X Axis Response 


0.08 


0.92 


Y Axis Response 


0.21 


0.77 


Z Axis Response 


0.21 


0.77 


elevator 


12.1 


1.0 


aileron 


12.3 


1.0 


throttle 


12.4 


1.0 


others 


0.62 


0.79 




1.90 


0.35 




2.17 


1.0 




2.35 


0.18 




2.26 


0.68 




4.54 


1.0 




5.86 


0.74 



elevator 




Frequency (rad/sec) 




Frequency (rad/sec) 



Figure 3.4: Control-Loop Bandwidth: Elevator Channel 



The response to two types of trajectory commands is of interest. The first 
is the reponse to a ramp co mman d in Y position. This corresponds to a change 
in the heading of the commanded trajectory. The response in terms of angle of 
bank and heading activity is shown in Figure 3.10. The controller achieves the 
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throttle 




Frequency (rad/sec) 




Frequency (rad/sec) 



Figure 3.5: Control-Loop Bandwidth: Throttle Channel 



aileron 





Figure 3.6: Control-Loop Bandwidth: Aileron Channel 

desired result of turning Bluebird to the required heading. The nonminimum 
phase response of the heading state is due to adverse yaw. 

The response to a ramp command in Z position corresponds to the re- 
sponse to a change in flight path angle of the commanded trajectory. Figure 3.11 
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X channel 




10' 2 1 0 1 10° io 1 io 2 



Frequency (rad/sec) 

Figure 3.7: Command-Loop Bandwidth: X Position Channel 



V channel 




Figure 3.8: Command-Loop Bandwidth: Y Position Channel 

shows the response to this command in terms of pitch angle, 6 , and altitude z. 
Note that the positive z direction is down, hence the negative pitch angle. 

Finally, the response to a constant wind disturbance is shown in Fig- 
ure 3.12. The wind orientation is such that it affects all three axes of Bluebird 
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Z channel 




-2 >1 O 1 2 

1 0 1 O 1 O 1 0 1 o 

Frequency (rad/sec) 



Figure 3.9: Command-Loop Bandwidth: Z Position Channel 
TABLE 3.5: COMMAND AND CONTROL BANDWIDTHS 



Loop 


Break Frequency 


Control Loop 


rad/sec 


Elevator 


4.0 


Aileron 


5.0 


Throttle 


1.0 


Command Loop 




X Axis 


0.75 


Y Axis 


1.0 


Z Axis 


1.0 



equally. The wind disturbance begins at t = 0. 
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Amplitude 




Figure 3.10: Bank Angle and Heading Response to Ramp in Y Command 




Figure 3.11: Pitch Angle and Altitude Response to Ramp in Z Command 
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Figure 3.12: Trajectory Error Due to a Constant Wind Disturbance 
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IV. CONTROLLER IMPLEMENTATION 
ON THE NONLINEAR PLANT 



In Chapter III, a linear controller was designed to control the trajectory 
of an air vehicle. However, this controller cannot be implemented, as is, on 
the nonlinear plant. It would only be effective for commanded trajectories that 
represent relatively small perturbations from the specific trajectory for which 
it was designed. Intuitively, it is clear that the dynamics of the plant do not 
depend on the heading angle 'P. Furthermore, the orientation of the force of 
gravity is the only change in the dynamics of the air vehicle due to changes in 
$ or 0. It turns out that these issues were addressed in [Ref. 10], where a new 
methodology for implementing controllers on nonlinear plants is proposed. The 
method involves differentiating some of the inputs to the controller, hence the 
term, ©-implementation. 

This chapter begins with a general description of the structure of ©- 
implementation. Furthermore, the specifics of its implementation on the non- 
linear model in SIMULINK are discussed. Next, the fidelity of the nonlinear 
simulation is improved by incorporating output feedback. This step involves 
inclusion of high fidelty sensor models and Kalman filters. Finally, all of the 
pieces of the complete nonlinear simulation are brought together in SIMULINK. 

A. ©-IMPLEMENTATION 

Using the development in ^Chapter II, the vehicle dynamics can be ex- 
pressed in state space form as follows: 
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— B v cg = f v ( B v cg B uj B ,u) +f R{A)g 
j t b ub = M B v cg , B 
j t P, * =' ft(A)% 

^ A = 5(A) b w b , 



(4.1) 



where f v and f w axe continuously differentiable and u € ft 6 denotes the vector 
of control inputs. To condense the notation, we define 



x„ := 



u eg 

3 ub 



x p := 



1 eg 

A 



/*(•) == 



?K)9 

0 



L(.) := 



/i(0 := 

K) 



/.(•) 

/«(•) 



0 



0 

5(.) 



where xi € ft 6 , x 2 € ft 6 and L 6 .ft 6 * 6 . Furthermore, we define 

€ ft 6 , 



r := 



ftc 

A c 



(4.2) 



Q :={ 



as the vector of linear and angular position commands that Bluebird must 
track. With this notation, the dynamics of the augmented plant can be written 
as follows: 

' = /i(x„,u) + / 2 (x p ) 

Xp L(x p )x v 

e i = [I 0](y 2 - r) ( . 

e 2 = [0 r\{y 2 - r) 1 ' 

Vi = h(x v ,u) 

V2 = Xp 

where y\ and y 2 axe the available measurements, e\ and e 2 are the trajectory 
errors separated into linear and angular components. Notice that L is only a 
function of the orientation vector A = [0 I]x p . For simplicity of exposition, we 
have not included any extra dynamics for the actuators or sensors. 
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The set £ of trajectories where the plant Q is expected to operate is given 



by 



£ • — { (^vo > ®po 5 r o) • 



X t>t) — 



Vo 

u; 0 



Xpo = r o 



> r o € r), 



. fl{ x voi x poi u o) "l" fi{ x po) — 0 . 

where T corresponds to the set of prescribed linear and angular positions of 



Bluebird. This is a broad definition of trimmed flight. Notice that it does 
not preclude the presence of an inertial acceleration due to centripetal force. 
As usual, we restrict the angular positions to some subset of [— 7r/2, tt/2] x 
(— 7t/2, 7t/2) x [— 7 t/ 2, 7r/2], as the inverse of L(.) is not defined at 9 = ?r/2. 
Notice, from the definition of x ^ € £ and equation (4.1) it follows that: 



B Wfl G £ —* B ub = constant 
A G £ — > A = constant. 



Notice that the set £ is easily parameterized by = r 0 € I\ Given 
u 0 , r 0 ) G £, we obtain 



yio '■= 

3 / 2 0 == x po 

e io : = [I 0](y 2o — r o) := 0 
e 2 0 := [0 /](y2 0 - ri)) := 0.. 

Let 8x v , 8x p , 8u, 8yi, 8y 2 , 8r, 8e 1 and 8e 2 correspond to small per- 
turbations of x„, x p , u, j/x, j/2, r, ex, and e 2 about the nominal values 
Xu,, , Xpo , ito, yi 0 , t/ 2 o, r 0 , ex 0 , and e 2o respectively. The family of lin- 
earized models associated with the rigid body Q and the set £ is defined as 
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Qi := {Qi (r 0 ), r 0 € T}, where 



Qi ( r o) := - 



Sly 

6 x v 

8yl 

hi 

8c\ 

8t 2 



«0 



and 



A\8x v + A 2 8x v + B8u 
L8x v + A^Sip 
C\8x v -f- D\8\i 
8 x p 

[/ 0](i!/ 2 - Sr) 

(0 /](«!), -Sr) 
v(i w ,r 0 ,u 0 ), 



(4.4) 



M 



0 gxf i2(Ao)5- 1 (Ao) 

0 0 



A,= 



0 — gi2(A 0 )Vo x 5 -1 (A 0 ) 

0 0 

(4.5) 



Matrices Ai and A 4 were derived using the identity in Appendix C. The intent of 
this derivation is to isolate the plant dynamics that axe a function of Ao. Notice 
that A-i represents the contribution of the force of gravity to the dynamics of 
Bluebird and A\ represents the sensitivity of Bluebird’s trajectory to changes 
in the air vehicle’s spatial orientation. 



Let r 0 € S be given. Define 



Ok ( r °) = < 



8 x v = Ai 8 x v + A 2 8 x p 4- B 8 u 
8 x p 8 xy 4 A^ 8 xp 
C\ 8 x v + D\ 8 u 



hi 

hi 

8c 1 
8c 2 
a 0 



= 8ip 

= [I 0](<fy 2 - 8 r) 
= [0 I]( 8 y 2 - 8 r) 

= u(x^,r 0 ,u 0 ), 



(4.6) 



where 



A 2 



0 g xf R( A 0 ) 

0 0 



A4 = 



0 — vox 
0 0 



(4.7) 



Notice that Equation 4.6 decribes the linear model of Bluebird used for 
the design of the LQR controller in the previous chapter, where ro was chosen 
as the origin of {/} with A 0 equal to zero. Note, at this condition, f i2(0) = I. 
Recall, the structure of the controller developed in Chapter III is given by: 
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ICt = 



(4.8) 



{ <5*cl — B c2 8l c2 

Sx c2 = [Se x 0] T 

Su = Cc\Sx c \ + C&Sxc 2 A D c \8yx + [D& D c2 \[8e\ 6t 2 ] T . 

Based on fCi , we propose to implement the following controller for the 
nonlinear plant Q : 



K (A) := < 



id = B c2 L~\ A)[ ei 0] r 

x<a = CciXci + C c2 L~ l {K)[ex 0] T + D cl yx 

+D c2 L-'( A)[0 e 2 ] T 
A = [0 1]y 2 

« = x c2 + D c3 L~ 1 (A)[ex 0] T . 



(4.9) 



Figure 4.1 shows the block diagram of fC (A). Notice that A serves as a 
gain scheduling variable. As a function of A, L~ l serves to properly resolve 
the trajectory error at the input to the controller. Note, the controller forms 
the derivative of the measured outputs, y\. Recall, y\ is the measurement of 
the states x v and the dynamics of x v axe essentially independent of the air 
vehicle’s spatial orientation or position in {/}. An integrator at the output of 
the controller serves to recover the properties of the linear design. The error is 
formed using the outputs y 2 and the commanded trajectory r. Recall, y 2 is the 
measurement of the states x p . L' 1 serves to resolve this error, originally formed 
in {/}, in {5}. 

It turns out that the implementation of Figure 4.1 has an important prop- 
erty discussed next. First we need to make the following assumptions: 

A\. Dim(x c2 ) — dim(u) — dim(y 2 ). 

A2. The matrix 



si - A c x B c2 

—Crt C c 2 



has full rank at s = 0. 



45 




0 


0 Bc2 


C c i C c 2 


Del [-Dc3 Dc2] 



Figure 4.1: Block diagram of the nonlinear controller K. (A) 

A3. The matrix pair (A e i,C c i) is observable. 

Also denote by T(Qi 0 tCi) : 6r — ► Sy? the closed loop linear system that 
results from connecting Qi 0 to K\ , and by T(Qi 0 tC{ )(s) its transfer matrix. 
Similarly, we let T\{Q , K )(tq) denote the linearization of the closed loop system 
T(Q , tC) at the equilibrium point determined by ro and let Ti(Q , K, )(ro)(s) 
be its transfer matrix. Then the following hold. 

• the feedback systems 7 }(Q , K, )(ro) and T(Qi 0 , fCi ) have the same closed 
loop eigenvalues; 

T,(5 , K )(r 0 )(s) = L(A 0 )T(g lt , K, )(s)i-‘ (A 0 ), 

Ao = [0 /]z M 

Thus, the eigenvalues of the linearizations at each operating point are 
preserved; furthermore, the input-output behavior of the linearized operators 
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is preserved in a well-defined sense. The reader is referred to [Ref. 10] for a 
complete discussion on approximations to this method that avoid using pure 
differentiation. The proof of this result is contained in Appendix C. 

B. D-IMPLEMENTATION OF THE CONTROLLER IN SIMULINK 

In general, three steps are required to implement the linear controller 
developed Chapter III on the nonlinear plant. First, the controller needs to 

form the derivative of x v as per Equation 4.9. Recall, x v is equal to the vec- 

• • • 

tor [ B Vcg, B u>B]‘ The first three elements of the vector, B v cg , axe available as 

processed acceleration outputs from the IMU. Therefore, the controller need 

not compute the derivative of B Vcg since it will have it available directly. The 

derivative of b u>b is computed by the controller. 

Secondly, the controller needs to act on the vectors, t\ and t-i. The linear 
position error is formed as the difference in commanded and current posi- 
tion. The vector t\ is then multiplied by the transformation matrix B R. This 
effectively resolves the linear trajectory error in the body-fixed reference frame. 
Along this position trajectory, there exists a corresponding trajectory of Euler 
angles. Recall that x v = constant is one of the constraints on the set of tra- 
jectories. This includes a broad range of flight conditions such as steady turns, 
steady pull-ups, climbing or descending turns, or constant heading. While it 
is natural to define the linear trajectory as a sequence of positions, it is more 
convenient to define the derivatives of the Euler angles rather than the values of 
the angles directly. Consider, for example, that for many trajectories of interest 
in £ , the components of A can be described by: 6 = 0; 0 = 0; and 'P = desired 
turn rate. Furthermore, the relationship between the rates of change of the 
Euler angles and b ub is used as a means of resolving the angular position error 
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in {5}. Recall, 



b u b = S~\ A)A. (4.10) 

Therefore, the derivative of the Euler angle states is formed and the commanded 
Euler angle rate is removed to form e' 2 . This error is resolved in {5} using 
Equation 4.10. The integrator at the end of the controller recovers the effect of 
the initial differentiation. 

Thirdly, the required error states axe formed by integrating the rotated 
linear trajectory error vector ej. Figure 4.1 indicates that integral action is 
accomplished at the output of the controller in order to recover the original 
properties of the linear design. This accounts for one of the integration steps 
on the error. Therefore, only one additional integrator is required to provide 
double integration action on the trajectory error ej. 

Figure 4.2 shows the P-implemented controller in SIMULINK, file plantl.m . 
See Appendix B for a complete description. The LQR gain has been parsed into 
several separate matrices for clarity of control action. 

C. GENERATION OF THE TRAJECTORY COMMANDS 

The commanded trajectory is specified with respect to the inertial refer- 
ence frame. At this point, it is assumed that a knowledge of the air vehicle’s 
performance capabilities is known and that the specified trajectory is within 
those capabilities provided there is no wind. The air vehicle has certain airspeed 
restrictions with respect to the air mass that cannot be violated regardless of 
the desired trajectory to be tracked. These restrictions typically provide lower 
and upper bounds on the velocity of the air vehicle with respect to the airmass. 
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Figure 4.2: ©-Implementation of controller on Bluebird 

An example of a lower limit is the stall speed of the air vehicle. Such limits are 
usually based on fundamental physical limitations of the airframe and a "fly- 
able” trajectory can become ”unflyable” under certain conditions. Therefore, 
a logic block positioned between the commanded trajectory and the controller 
ensures that the commanded trajectory can be flown at current flight conditions 
within user defined indicated airspeed limits, shown in Figure 4.3. The com- 
manded linear trajectory enters the block as a time stamped position fix in the 
inertial reference frame. Onboard sensors provide both inertial velocities from 
the IMU and air mass velocities from the pitot-static system. Furthermore, a 
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dual vane device instrumented on Bluebird provides both a and /3 measure- 
ments. Note: a close approximation to these readings could be obtained from 
IMU measurements as outlined in Chapter II. 




“K? 



Figure 4.3: Commanded Trajectory Logic Block 



With these measurements, the wind vector resolved in {/} is calculated as: 



'W =’ B R b v v ~' b R$,R w v c , 



(4.11) 



where 



' v t ■ 

= 0 

. 0 . 

and V t is the indicated airspeed obtained from the pitot-static system. 

The commanded indicated airspeed of the air vehicle, Vt-cmd , is calculated as: 
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where / v cm( i is the numeric derivative of the commanded linear trajectory. If the 
commanded indicated airspeed is not within specified limits, the commanded 
trajectory is altered as follows. The angles 6 and rp are calculated as follows: 



6 = tan 1 (v T /v x ) (4-12) 

and 

rp = sm~ 1 (u v /| / u cmt/ |) (4.13) 

where the components of I v cm d are [u x , v y , v z ] T . Note that the angles 9 and rp 
define the commanded velocity vector’s orientation in {/}. 

Finally, the amount that Vt-and is outside of indicated airspeed limits is 
subtracted from the magnitude of I v cm d, resulting in the magnitude of the new 
commanded inertial velocity, termed V. V is then resolved in {/} according to: 



Vmod — 



cos(6)cos(tp) —cos(6)sin(ip) —sin(9) 
sin(ip) cos(xp) 0 

_ sin(9)cos(rp) —sin(9)sin(rp) cos(9) 



(4.14) 



where ! v mo d is the new commanded inertial velocity. This command is inte- 
grated and sent to the controller as the commanded trajectory. The MATLAB 
.m file that implements this logic is windlogic.m and can be found in Appendix 

A. 



The net effect of the trajectory logic block is simple. When Bluebird runs 
up against one of its airspeed limits, the commanded trajectory can no longer be 
followed. A choice is made to maintain the direction of the commanded velocity 
but change its magnitude. Notice that this method does not affect the turn rate 
associated with the trajectory and subsequently, no processing of the angular 
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Bluebird Heading - degrees 



trajectory commands is required. In this way, Bluebird is never commanded to 
fly a trajectory that would force it to exceed the performance limits. 

The performance of the trajectory logic block can be seen in Figure 4.4. 
The lower limit for Bluebird was arbitrarily choosen as 63 feet per second. Blue- 
bird is initially flying due north at a ground speed of 73 feet per second, crabbed 
into 20 feet per second of wind from the west. The commanded trajectory turns 
90 degrees to the east. Notice that the original trajectory would result in an 
indicated airspeed of 53 fps while the revised trajectory results in a commanded 
trajectory of 63 fps. 





Figure 4.4: Example of Commanded Trajectory Revision 

If the commanded trajectory is generated using a velocity rather than 
position schedule, the differentiation block in Figure 4.4 can be removed. A 
velocity schedule is specified in {/} as a sine wave of appropriate magnitude, 
frequency, and phase along the x and y axis. The commanded ground speed 
corresponds to the magnitude and the commanded turn rate corresponds to 
the frequency of the sine function. Constant heading flight is a subset of these 
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trajectories where the frequency is zero. This method of generating trajectory 
commands makes determining turning rate (A and) simple. As an example, con- 
sider the case of generating the velocity schedule for a circular flight pattern 
at a ground speed of 100 fps and a turn rate of 0.1 radians per second. The 
derivative of the commanded trajectory, r, parameterized by time is: 



x axis velocity 
y axis velocity 
z axis velocity 
$ 

0 



100sm(0.1t + pi/2 ) 

100sm(0.1t + 0) 
desired climb or descent rate 
0 
0 

0.1 



The SIMULINK block that generates the commanded trajectory is shown 
in Figure 4.5. 




upper 

limit 



Figure 4.5: Generation of Trajectory Commands 
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D. STATE FEEDBACK TO OUTPUT FEEDBACK 



Up until this point, the development of the tracking controller has dealt 
exclusively with full state feedback. This section will detail a method whereby 
the full state feedback is replaced by high fidelity models of the onboard sensors 
in conjunction with Kalman filters designed to provide optimal state estimates, 
using onboard sensor data. 

1. Sensor Modeling 

This work builds upon sensor models developed in two prior theses. 
In [Ref. 1], a detailed model of the inertial measurement unit, IMU, is devel- 
oped. In [Ref. 2], a detailed model of the GPS unit is developed. The IMU is 
a compact, lightweight, low power unit which integrates nine sensor measure- 
ments in one box. These sensors are three axis accelerometers, three axis rate 
gyros, pitch and roll inclinometers, and a magnetometer. The accelerometers 
are instrument grade, signal conditioned, and temperature compenstated. Full 
scale output is + / — 3 g’s. The accelerometer’s frequency response is flat past 
100 Hz. However, the antialiasing inside the IMU limits the effective bandwidth 
of all of the sensors to 20 Hz. An internal initialization program allows the unit 
to compenstate for accelerometer bias and cross axis error. Table 4.1 shows the 
specifications of the accelerometers incorporated into the sensor model [Ref. 8]. 

The rate gyros used by the IMU are solid state vibrating element 
angular rate sensors. This relatively new technology uses no moving parts. A 
piezoelectric bender element is mounted end to end but rotated at a 90 degree 
angle. The element fastened to the base is resonantly driven such that the 
sense element swings a reciprocating arc. Under zero angular rate conditions, 
the motion of the sense element due to the drive element does not produce any 
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TABLE 4.1: ACCELEROMETER CHARACTERISTICS 



Acceleration Range 


±2g's 


Acceleration Bandwidth 


20 Hz 


Acceleration Bias 


0.2% of Full Scale 


Acceleration Scale Factor 


0.2% of Full Scale 


Acceleration Noise Floor 


0.0005 g's 


Cross Axis Sensitivity 


0.5% of Full Scale 



bending of the sense element. When a rate of rotation exists, Coriolis forces 
cause momentum to be transfered into the plane perpendicular to the motion 
of the drive element, thus causing bending of the sense element. A pressure 
transducer picks up a signal from the sense element when it is bent that is 
proportional to the angular rate with a phase dependent on the direction of 
rotation. Figure 4.6 shows the configuration of two rate sensors mounted in a 
"tuning fork” configuration. Table 4.2 shows specifications of the rate sensors 
incorporated into the sensor model [Ref. 8]. 




Figure 4.6: Angular Rate Sensor 




TABLE 4.2: ANGULAR RATE SENSOR CHARACTERISTICS 



Rotational Rate Range 


±114.6de<7/sec 


Rotational Rate Bandwidth 


20 Hz 


Rotational Rate Bias 


2.0% of Full Scale 


Rotational Rate Scale Factor 


0.5% of Full Scale 


Rotational Rate Noise Floor 


0.05% of Full Scale 


Cross Axis Sensitivity 


0.5% of Full Scale 



The inclinometers utilize a .liquid crystal pendulous sensor. It is a 
low bandwidth sensor ( approximately 0.12 radians per second) that is meant 
to be integrated with the rate sensors for high bandwidth measurements of 
angular position. The fluxgate magnetometer provides heading measurements. 
The specifications incorporated into the Euler angle sensor models are shown 
in Table 4.3 [Ref. 8]. 

TABLE 4.3: INCLINOMETER AND MAGNETOMETER CHARAC- 
TERISTICS 



Pitch and Roll Range 


±50 deg 


Pitch and Roll Bandwidth 


1/2 Hz 


Pitch and Roll Accuracy 


0.2 deg 


Heading Range 


±180 deg 


Heading Accuracy 


3.0 deg 


Heading Repeatability 


0.5 deg 


Heading Linearity 


0.5% 



GPS provides data in a form that can be converted to local tangent 
plane position. A brief summary of the errors included in the GPS model 
follows. 

GPS Error Sources 
• Selective Availability 
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- intentional degradation of pseudorange signal by Departmant of Defense 

• Clock Differences 

- drift and bias in GPS clock 

• Ephemeris Error 

- error introduced in converting pseudoranges to inertial position fix 

Each of these sensor components is simulated in block diagram form in SIMULINK 
utilizing internal modeling principles based on manufacturer specifications and 
known sources of error. The upper level SIMULINK diagram of these sensor 
models is shown in Figure 4.7 and contained in the SIMULINK file plantl.m. 




Figure 4.7: Sensor Models in SIMULINK 
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2. Kalman Filtering 



With the outputs of the modeled sensors available, Kalman filters 
were developed along the linear and angular position channels to provide optimal 
estimates of the states for the controller. An approach analogous to the LQR 
design was used. Consider the linear system described by: 

x = Ax + Bu + Gw /a ic\ 

z = Cx + v (4 - 16) 

where v and w axe zero mean white noise with respective power spectral densities 
of V and W. 

A gain matrix, L, was found such that the Kalman filter given by: 



x = Ax + Bu + L(z — Cx) (4-17) 

produces an optimal estimate of x. The Kalman gain L is calculated as follows: 

L = YC t V~\ 

where Y is positive semidefinite and solves the algebraic Riccati equation: 

AY + YA t - YC t V~ x CY + GWCf = 0 

A synthesis model was formed that included the dynamics of the 
original plant. The IMU used in Bluebird incorporates an initialization routine 
that removes steady state bias from the sensors. Therefore, extra dynamics were 
not required in the Kalman filter to compensate for bias. The design process was 
primarily driven by the bandwidth limitations of the inclinometers and GPS. 
The values of V and W were used as ” knobs” in the iterative design process. 
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For the GPS sensor, a break frequency of 1/2 radian per second was desired. 
For the Euler angle sensors, a break frequency of 1/10 radian per second was 
desired. It was also desired to wash out the accelerometers and rate sensor 
biases at low frequencies. 

The Kalman filter for the position estimate blends processed ac- 
celerometer outputs with the GPS position fixes. Figure 4.8 shows the frequency 
response along the x axis of the Kalman position filter. The other two axes have 
similar dynamics. 




Figure 4.8: Frequency Response of Position Filter 

The Kalman filter for the Euler angles blends the outputs of the angu- 
lar rate sensors, converted to Euler angle rates, with Euler angle measurements 
from the inclinometers and magnetometer. Figure 4.9 shows the frequency re- 
soonse of one channel of the Euler angle filter. The remaining two angle channels 
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axe similar. 




Figure 4.9: Frequency Response of Euler Angle Filter 

E. INTEGRATION OF THE FULL NONLINEAR SIMULATION 

The full nonlinear simulation can now be pieced together. Recall in Chap- 
ter II, the nonlinear rigid body dynamics were implemented in a SIMULINK 
block labeled Equations of Motion. If the simulation is expanded to include the 
effects of a moving airmass, the dynamics of Bluebird can be simulated at an 
arbitrary flight condition. Wind is usually referenced with respect to the iner- 
tial reference frame, therefore a SIMULINK block, Wind , is included in the full 
simulation whose output is a vector w v comprised of the wind velocity resolved 
in {/}. Next, the wind velocity is resolved in {13} and added to the inertial 
velocity of Bluebird, B v cg . The result is the velocity of Bluebird with respect 
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to the airmass, resolved in {B}. This velocity, rather than B v cg , is used when 
computing the aerodynamic contribution to the total external force and moment 
used in Equation 2.25. A more detailed description of how this is accomplished 
in the .m file state.deriv.m is contained in Appendix A. 

All sixteen states axe sent to a SIMULINK block, Sensors, that models 
the avionic sensor suite onboard Bluebird. The output from these sensors is 
appropriately processed in a SIMULINK block, Kalman filters. The filtered 
output is directed to the T> -Implemented Controller block. The commands to 
the controller come from a trajectory block which uses the measured outputs 
from the filter block to process the commanded trajectory. The controller gen- 
erates actuator commands necessary to maintain the air vehicle on the com- 
manded trajectory, thus completing the loop. The complete top level view of 
the SIMULINK nonlinear simulation is shown in Figure 4.10 and contained in 
the SIMULINK file plantl.m. 
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Figure 4.10: SIMULINK Diagram of the Full Nonlinear Simulation 
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V. APPLICATION TO THE CONTROL OF 

BLUEBIRD 



A. ©-IMPLEMENTED CONTROLLER PERFORMANCE CHAR- 
ACTERISTICS 

The performance characteristics of the trajectory tracking controller de- 
signed in Chapters III and IV were evaluated using a two step process. First, the 
sensor and filter blocks were removed and the controller was connected directly 
to the nonlinear equations of motion block. Utilizing pure state feedback, Blue- 
bird was flown along two fundamentally different types of trajectories. These 
two trajectories served as general examples of the set of all trajectories defined 
in Chapter IV, Equation 4.4. Next, the sensor and filter blocks were added. 
A general deparure and arrival trajectory, which is a combination of the two 
types of trajectories tested in the first step, was commanded and flown with 
the controller utilizing output feedback. Data from this simulation were used 
as input to a virtual prototype simulation discussed later. 

The dynamic flight simulations were started using the same initial condi- 
tion. At this initial condition, Bluebird is aligned with the positive x-axis and 
trimmed for level flight at 73 fps. The positive x direction is considered to be 
heading north. The mechanics of the dynamic simulation use a right-hand or- 
thogonal coordinate system described in Chapter II. As such, the positive y-axis 
is pointing east and the positive z-axis is pointing down. This choice is con- 
venient from a computational standpoint since it coincides with the body-fixed 
reference frame of Bluebird at the initial condition specified above. Typically, 
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however, the positive z axis is considered to be pointing up. If the right-hand 
orthogonality is maintained, the positive y-axis would then point towards the 
west. For ease of visualization, this is the cooridinate system used to display 
the results of the simulations. 

The simulations axe intended to evaluate the capabilities of the controller 
in terms of the nature of trajectories that can be followed in a stationary air mass 
and in an air mass moving at constant velocity. Two basic kinds of trimmed 
flight serve as the bases for the test trajectories. Each test trajectory is flown in 
no-wind conditions, and then again with the wind added at some point during 
the flight. 

The simplest form of trimmed flight is constant velocity, constant heading. 
This corresponds to a trajectory defined by a ramp command in inertial posi- 
tion. This was the basic trajectory that the controller was designed to track. 
Figure 5.1 shows a three dimensional plot of the first test trajectory. In this 
case the trajectory encompasses 30 seconds of flight heading north at 73 fps 
followed by a 90 degree turn to join a trajectory heading east at 73 fps while 
climbing at 300 feet per minute. On the second flight, a wind of 20 feet per 
second from the north is added at the time the turn is commanded (elapsed 
time = 30 seconds) . 

Figure 5.2 contains the first four graphs of flight data. The first graph 
shows the time history of Bluebird’s distance from the commanded trajectory. 
The next three graphs show the time history of the Euler angles. Consider the 
baseline flight (no- wind data). Bluebird begins the turn at an elapsed time of 
30 seconds and exits the turn at an elapsed time of 45 seconds. Approximately 
30 seconds later the trajectory error is nearly zero. In the presence of 20 feet 
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per second of wind from the north, the trajectory error also goes to zero in 
about the same period of time. The graphs of the Euler angles indicate that 
Bluebird is flying wings level, crabbed into the crosswind, which is the desired 
result. Figure 5.3 shows the groundspeed and indicated airspeed during the 
flights. Notice that in both cases, the commanded groundspeed of 73 feet per 
second is eventually maintained. Finally, Figure 5.4 shows the time history of 
the control activity during the flights. 

Trimmed flight does not necessarily have b ub equal to zero. For instance, 
any steady turning manuever fits the definition of trimmed flight given in Chap- 
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Figure 5.2: Trajectory #1: Position Error and Euler Angles 



ter IV. Figure 5.5 shows the three dimensional plot of the second test trajectory. 
In this case the test trajectory is a helix flown at 73 feet per second. The turn 
rate is one revolution per minute and the helix angle is 4 degrees which corre- 
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Qroundspeed 





Figure 5.3: Trajectory #1: Velocity Data 

sponds to a climb rate for Bluebird of 300 feet per minute. For this test, no 
indicated airspeed limits were placed on Bluebird. 

Consider Figure 5.6 which shows the position error and Euler angle time 
history for the helix trajectory. Notice that with no wind the controller manuev- 
ers Bluebird to join the commanded trajectory with zero error in steady state. 
The constant pitch and bank angles confirm the steady state performance. On 
the second flight, a wind of 20 feet per second from the east was added at the 
start of the helical trajectory (elapsed time equal to 40 seconds). Intuitively, 
it is clear that the bank and pitch angle must vary as Bluebird traverses the 
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Figure 5.4: Trajectory #1: Control Activity 



helix. To an observer on Bluebird, the wind, while constant in {/}, represents 
a sinusoidal disturbance. This explains the sinusoidal nature of the position er- 
ror around the helix. Figure 5.7 shows the groundspeed and indicated airspeed 
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around the helix during the two flights. For the no- wind flight, the commanded 
groundspeed is maintained while for the flight in wind, a one foot per second 
oscillation is experienced. Finally, Figure 5.8 shows the time history of the 
control activity during the flights. 

Four additional flights were flown using the helix trajectory in an attempt 
to ascertain the sensitivity of the position error to changes in commanded turn 
rate and wind velocity. Figure 5.9 shows the position error around the helix 
for three different turn rates with a constant wind of 10 feet per second. The 
dashed line corresponds to a turn rate of 3 degrees per second or 2 minutes 
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Figure 5.6: Trajectory # 2 : Position Error and Euler Angles 



per revolution; the solid line corresponds to a turn rate of 6 degrees per second 
or 1 minute per revolution; and the dash dot line corresponds to a turn rate 
of 9 degrees per second or 40 seconds per revolution. The error increases with 
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Figure 5.7: Trajectory #2: Velocity Data 



increasing turn rate. The 9 degree per second turn rate corresponds to a steady 
state angle of bank of thirty degrees, no wind. It may not be desireable to 
co mm and trajectories requiring more than a certain angle of bank and this may 
place an upper bound on the error. 

Figure 5.10 shows the position error around the helix for three different 
wind velocities at a constant turn rate of 6 degrees per second. The wind varies 
in velocity from 10 to 25 feet per second. The local maxima values of the 
position error are proportional to the magnitude of the disturbance. 
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Figure 5.8: Trajectory #2: Control Activity 

B. AN AIRPORT DEPARTURE AND ARRIVAL FLIGHT SIM- 
ULATION 

In many cases, the trimmedflight condition of an air vehicle changes often. 
A good example of this would be a standard instrument departure or arrival 
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Figure 5.9: Trajectory #2: Position Error for Varying Turn Rates 



to an airfield. These trajectories are typically a combination of constant radius 
turns and wings level flight, while often climbing and descending. Trajectory 
position errors become critical when the air vehicle is on final approach with a 
constant heading, constant velocity trajectory. 

Consider Figure 5.11. If an airfield is imagined to be located at the origin, 
then this trajectory would be indicative of a typical departure followed by a 
typical arrival to that airfield. The scenario utilizes turning trajectories of three 
different radii connected by straight line trajectories. The commanded velocity 
is a constant 73 feet per second throughout. Wind is initially zero. Thirty 
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Figure 5.10: Trajectory #2: Position Error for Varying Wind Velocities 

seconds into the flight, the wind is added at 10 feet per second from the east. 
At 90 seconds into the flight the wind is increased to 20 feet per second from the 
east. Finally, with Bluebird on final approach tracking a 4 degree glideslope, 
the wind is rapidly shifted 90 degrees to the north and decreased in magnitude 
to 5 feet per second. 

Figure 5.12 shows the time history of the position error, wind velocity, and 
Euler angles during the flight. Figure 5.13 shows the time history of the control 
activity. Note, however, the relative difficulty of analyzing data of this nature 
as it is somewhat difficult to visualize. Flight simulation data was saved to a file 
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and processed for compatiblity as an input file for a 3-D visualiztion software 
package, Designer’s Workbench. A virtual prototype of Monterey Airport and 
Bluebird was developed in [Ref. 13]. The simulation was then run as a departure 
and arrival to the virtual prototype airfield. In Designer’s Workbench, the flight 
can be viewed from multiple perspectives and virtual prototypes of standard 
cockpit displays further enhance visualization. One possible result is captured 
on video tape, [Ref. 14]. 
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Figure 5.12: Airfield Scenario: Position Error and Euler Angles 
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VI. CONCLUSIONS AND 
RECOMMENDATIONS 



A. CONCLUSIONS 

Based on the data presented in this thesis, the following conclusions axe 
drawn. 

• SIMULINK provides an effective environment for the developement of non- 
linear simulations for air vehicles. As a result of this development, a linear 
model of the plant at an arbitrary trim condition is easily obtained for 
design purposes. 

• LQR design techniques utilizing a synthesis model and weighting ” knobs” 
provide a straight forwaxd means of obtaining satisfactory controller gains 
for MIMO systems while meeting design requirements. 

• ^-Implementation of the linear trajectory tracking controller allows the 
controller to operate effectively on the nonlinear plant. In no- wind flight 
conditions, trajectories defined by an arbitrary [u 0 , u>o] are tracked perfectly 
in steady state. For flight conditions with wind, rejection of a constant 
wind disturbance is accomplished along the family of trajectories defined 
by an arbitrary [uo,u>o = 0]. However, for turning flight, a constant wind 
disturbance in {/} is seen as a sinusoidal disturbance in {!?} and a sinu- 
soidal tracking error results. For moderate bank angles and turn rates, the 
errors axe usually small. 
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• Preprocessing of the trajectory commands by an adaptive filter allows for 
steady state control of the air vehicle’s velocity with respect to the air mass 
in the forward path, thus not affecting stability. With sufficient margins for 
transient deviations in indicated airspeed, this would ailieviate the major 
concern of stalling the air vehicle when tracking an inertial trajectory in 
wind. 

• When analyzing a nonlinear plant and controller, test simulations are vital 
and in some cases the only means of performance evaluation. The three 
dimensional plots and time history graphs are fine for simple trajectories, 
but are difficult to analyze for more complex cases. The capabilities of a 
virtual prototyping software package, like Designer’s Workbench, are im- 
pressive. The enhanced situational awareness and visualization capabilities 
of watching the designed controller operate on a virtual prototype allow 
for a "pilot’s perspective” feedback not otherwise attainable on the desk 
top. 

B. RECOMMENDATIONS 

Based on the conclusions presented above and the experience of developing 
the simulation package presented in this thesis, the following recommendations 
are made. 

• While the rigid body equations of motion are nonlinear with respect to 
the kinematics involved, they are completely linear with respect to the 
stability and control derivatives. The constant coefficient stability and 
control derivatives could be replaced by functions when further flight data 
is available. 
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• A similar design process used for the trajectory controller could be done 
using Woo design methodology. 

• The trajectory preprocessor could be used to convert an inertial fixed tra- 
jectory into an air mass fixed trajectory. This might have applications 
where the air vehicle’s inertial position is of secondary importance com- 
pared to its performance with respect to the air mass. 

• Running simulations real time in Desiner’s Workbench rather than using 
batch post processed data would be the next logical step. Further work 
might lead to virtual prototype visualizations based on real-time simula- 
tions or downlinks from actual air vehicles. 
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APPENDIX A: MATLAB FILES 



The SIMULINK models of Bluebird (plant l.m & plant2.m) use the fol- 
lowing MATLAB .m files as user defined functions. 

STATE JDERIV.M 



•/•/•/ * 1 * 1*1 



X 

X 

X 

X 



X 

X 



X 

X 



X 



X 



X 



X 

X 

X 

X 

X 



X 



Function to calculate derivative of u,v,w,p,q,r 
based on 

1: kinematics 
2: gravity 
3: propulsion 
4: aerodynamics 



X 

X 

X 

X 

X 

X 

X 

X 



Variables brought from workspace: */. 

X 

x - [contrl inputs, state variables(l - 9), wind vel] */, 
= (da,de,dr,dtrt ,u,v,w,p q,r, phi, theta psi, wind xyz)'/, 

y. 

Variables called from function "blue_data" X 



rho * air density 
b = wing span 
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X 

X 

X 



y. 

y. 



x 

x 

y. 

x 

y. 

x 

x 

x 

x 



c = wing cord X 

s = wing area X 

Cfo = Steady state force term X 

Cfu = Stability derivative for control inputs X 
m = airplane mass X 

lb ■ inertia tensor matrix (body frame) X 

To = Thrust scale term X 

Pe = Engine post ion matrix X 

X 

X 

X 



function accel = state.deriv(x) 



XXXXXX Function call to get the aircraft data 



[uO,wO,rho,Cfx,Cfo,Cfu,Cfxdot,s,b,c,m,Pe,To,Ib] = blue.data; 



XXXXXX separate the combined vector into separate elements 

u = [x(l) ; x(2) ; x(3)] ; 
dtrt » x(4) ; 

state = [x(5) ; x(6) ; x(7); x(8) ; x(9) ; x(10)] ; 
lambda = [x(ll); x(12); x(13)] ; 
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wind = [x(14); x(15); x(16)] ; 



W/tV/tW/' calculate velocity wrt the airmass and form state vector 
W/, , /,V/,V, that will be used to calculate the aerodynamic forces/moments 



ias = u + wind 

statel = [ias(l)-uO; ias(2); ias(3)-w0; x(8) ; x(9) ; x(10)] ; 
7.7.77.7.7. calculate total velocity, vt 
vt “ (ias(l)~2 + ias(2)~2 + ias(3)“2)“ .5; 

77.7777. calculate qbar 
qbar = .5*rho*(vt''2) ; 

7.7.7777 calculate Ml 

Ml = diag([l/vt, 1/vt, 1/vt, (.5*b)/vt, (.5*c)/vt, (.5*b)/vt]); 
mm calculate M2 

M2 = diag([0, 0, (.5*c)/(vt“2) , 0, 0, 0]); 
mm, calculate Sprime 
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Sprime = diag([-l, 1, -1, b, c, b]*s); 



1XV/.V/, calculate Mu 

Mu = inv( [eye(3)*m,zeros(3) ; zeros (3) ,Ib]) ; 

calculate Tu2b 

alpha = ias(3)/vt; 
beta = ias(2)/vt; 

T1 = [cos (alpha), 0, -sin (alpha) ; 0,1,0; sin(alpha), 0, cos (alpha)] 
T2 = [cos (beta), -sin(beta) , 0; sin(beta) , cos (beta), 0; 0,0,1]; 
Tw2b = [T1*T2, zeros (3); zeros (3), T1*T2] ; 



mm. calculate Chi 

Chi = eye(6) - Mu*Tv2b*qbar*Sprime*Cfxdot*M2 ; 

'ItWItVIX calculate Propulsion matrix 

Prop = [ eye (3); 

0,-Pe(3),Pe(2); 

Pe(3),0,-Pe(l); 

-Pe(2) ,Pe(l) ,0] ; 
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’/.'/.'/.'/.V.'/. calculate gravity vector and rotation matrix {1} to {B} 

Rot = [1, 0, -sin(lambda(2) ) ; 

0, cos(lambda(l) ) , cos(lambda(2))*sin(lambda(l)) ; 

0, -sinClambda(l)) , cos(lambda(2))*cos(lambda(l))] ; 

Ru2b = [Rot ; zeros (3)] ; 

g = [0; 0; 32.174]; 

FgU - m*g; 

'ItVIXVL put the components due to gravity; thrust; and control surface 
% IXVI*Vt* deflections together for their contribution to x-dot 

thrust = Prop*To*dtrt ; 
gravity = Ru2b*FgU; 

Ctrl = qbar*(Tw2b*(Sprime*(Cfo + (Cfu*u)))); 
xdotu = (Mu*(ctrl+thrust+gravity)) ; 
nUU calculate kinematic contribution 

omegax = [0, -state (6) .state (5) ; state (6) ,0, -state (4) ; -state (5) .state (4) ,0] ; 
vxlb = (-inv(Ib))*(omegax*Ib) ; 
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Rot = [-omegax, zeros (3); zeros (3), wxlb] ; 



xdotrot = Rot* st ate; 



mU state vector feedback component xdot 



xdotcfx = qbar*(Mu*(Tv2b*(Sprime*(Cfx*(Ml*statel))))) ; 



mm add three components of xdot together and premult by inv(Chi) 



xdot** inv(Chi)*(xdotrot+xdotcfx+xdotu) ; 



mm return xdot 



accel**xdot ; 



EULER-B2I.M 



X 



X 



*/, Transformation [p q r] to lambda-dot 



7 . 



7 . 



7 . 



xxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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function ldot = eul_b2u(x) 



WltltVIt seperate the composite vector 'x' into [p q r] 
VfXVh and [phi theta psi] . 

omega = [x(l) ; x(2) ; x(3)] ; 
phi=x(4) ; 
theta= x(5) ; 
psi=x(6) ; 



mU calculate the rotation matrix {1} to {B} 
Vlt'IX'l* based on euler angles 

Rb2u = [1, sin(phi)*tan(theta) , cos(phi)*tan(theta) ; 
0, cos(phi), -sin(phi); 

0, sin(phi)*sec(theta) , cos(phi)*sec(theta)] ; 

'ItWItVI* calculate lamda-dot 
ldot = Rb2u*omega; 



EULERJ2B.M 



• /• /• h /• /• /• /• /i /• /• /• /•/•/•/•/•/•/•/•/•/•/• /# /• /• /• /• /• /• /• /• /• 
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y. 

y. 

x 



Transformation lambda-dot to [p q r] 



% 

X 

y. 



function omegadot = eul_i2b(x) 

IUU separate tbe composite vector ’x’ into lambda-dot 
U%U and [phi theta psi] . 

ldot = [x(l); x(2) ; x(3)] ; 
phi=x(4) ; 
theta® x(5) ; 
psi*x(6) ; 

IUU calculate the rotation matrix {B} to {1} 

UU based on euler angles 

Rb2i = [1, sin(phi)*tan(theta) , cos (phi) *tan (theta) ; 

0, cos (phi), -sin(phi) ; 

0, sin (phi)* sec (theta) , co s (phi) *sec (theta) ] ; 



UU calculate lamda-dot 



omegadot = inv(Rb2i)*ldot ; 
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P0STTI0N-B2I.M 



/ •/ •/ •/ •/ •/ 



X 

X 

X 

X 

X 

X 

X 

X 

X 

X 



X 

From the workspace: X 

X 

1: free vector ’u' resolved in {B} e(l:3) X 

2: euler angle vector {phi, theta, psi} e(4:6) X 

X 

Returns : X 

X 

1: free vector 'u resolved in {1} X 

X 



axxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



function ans = pos_b2i(e) 

XXXXX this will rotate the trajectory error through phi, theta, psi 
XXXXX from {b} to {i}. (3-2-1 transformation 



phi-e(4) ; 
theta=e(5) ; 
psi»e(6) ; 
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m_psi = [cos(psi) ,sin(psi) ,0 
-sin(psi) ,cos(psi) ,0 

0,0,1]; 

m_theta = [cos (theta) ,0,- sin (theta) 

0,1,0 

sin(theta) ,0, cos (theta)] ; 

m_phi = [1,0,0 
0 , cos (phi) , sin (phi) 

0, -sin (phi) ,cos(phi)] ; 



rotb2i = inv(m_phi*m_theta*m_psi) ; 
u e [e(l) ; e(2) ; e(3)] ; 
ans = rotb2i*u; 



POSITION J2B.M 



m 



v m 



•/. 
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y. 

y. 

y. 

y. 

x 

x 

•/. 

y. 



From the workspace: X 

X 

1: free vector 'u' resolved in {1} e(l:3) X 

2: euler angle vector {phi, theta, psi> e(4:6) X 

X 

Returns : X 

X 

1: free vector 'u resolved in {B} X 



y. 



x 



• /•/•/•/•/•/•/• 



function ans = pos_i2b(e) 

XXXXX this will rotate through phi, theta, psi 
XXXXX from {i> to {b}. (3-2-1 transformation) 



phi=e(4) ; 
theta=e(5) ; 
psi=e(6) ; 

m_psi * [cos(psi) ,sin(psi) ,0 
-sin(psi) ,cos(psi) ,0 
0 , 0 , 1 ]; 

m.theta = [cos (theta) ,0,-sin(theta) 

0 , 1,0 
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sin (theta) ,0, cos (theta)] ; 



m_phi = [1,0,0 
0,cos(phi) ,sin(phi) 
0,-sin(phi) ,cos(phi)] ; 



roti2b « (m_phi*m_theta*m_psi) ; 
u * [e(l) ; e(2) ; e(3)] ; 
ans = roti2b*u; 



LIMIT-LQGIC.M 



•/ •/ •/ •/ •/ •/ •/ •/ *I*J*I*I *1 v •' 



•X 



•/. 



% funtion to limit trajectory commands, if required '/, 
'/. •/. 
y, from workspace: ’/. 
'/. 1: commanded inertial velocity */, 
'/, 2: inertial wind '/. 
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3: lower IAS limit 



X 

X 



X 

X 

% 

% 



4: upper IAS limit 



7 . 



returns: revised commanded velocity 



X 



y. 



x 






.xxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



function vcom=limit (u) 



y.7.7.7.7. seperate u 



vel.i = [u(l> ;u(2) ;u(3)] ; 
wind.i = [u(4) ,u(5) ,u(6)] ; 



ll-u(7) ; 
ul=u(8) ; 



XXXXX calculate magnitud and direction of commanded velocity 



gs=sqrt(vel_i(l)“2 + vel.i (2) ~2) ; 



ang=atan2 (vel.i (2) , vel.i (1)) ; 



XXXXXX calculate commanded IAS (steady state) 
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vt= sqrt((vel_i(l)+wind_i(l))‘2 +(vel_i(2)+wind_i(2))“2 +(vel_i(3)+wind. 



'/.y.y.'/.y.'/. Prepare return variable (may not be limited) 
vcom = vel.i; 

mm Check limits and revise if outside 

if vt > ul; 
over = vt - ul; 

vcom(l) - (gs - over)*cos(ang) ; 

vcom(2) = (gs - over)*sin(ang) ; 
end; 

if vt < 11; 
under = 11 - vt; 

vcom(l) = (gs + under) *cos(ang) ; 
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vcom(2) = (gs + under) *s in ( ang) ; 
end; 



BLUEBIRD -DATA.M 




7. 



7. 



Aircraft data for Bluebird 



7. 



7 



7. 




function [uO,wO,rho,Cfx,Cfo,Cfu,Cfxdot ,s,b,c,m,Pe,To,Ib] = blue.dat a 

7.7.7.77. trimmed flight speed and angle of attack 

uO = 73.3; 
wO ® 0; 

7777.77. Density: Sea level- std day 
rho = .0023769; 

77777 derivative matrix due to state variables 
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mU rows: [CD CY CL Cl Cm Cn] 

7,UU col: [u v/U w/U p q r] 



Cfx = [00 .188 0 0 0; 

0 -.31 000 .0973; 

0 0 4.22 0 3.94 0; 

0 -.0597 0 -.363 0 .1; 

0 0 -1.163 0 -11.77 0; 

0 .0487 0 -.0481 0 -.0452]; 

*/>UU derivative matrix due to control inputs 

XXXXX rows: [CD CY CL Cl Cm Cn] 

mu col: [elev rud ail] 

Cfu - [.065 0 0; 

0 .0697 0; 

.472 .0147 0; 

0 .0028 .265; 

-1.41 0 0; 

0 -.0329 -.0347]; 



derivative matrix due to x-dot (alpha.dot k Beta.dot) 
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Cfxdot = [000000; 

0 0 0 0 0 0 ; 

0 0 1.32 000; 

000000 ; 

0 0 0 0 0 0 ]; 

U.U'/« steady state force vector 

Cfo = [.03; 0; .385; 0; 0; 0]; 
mU. physical dim. 

mu. WT =55 LBS. 

m * 1.7095; 
s = 22.38; 
b = 12.42; 
c = 1.802; 

UUX engine data (4 HP motor) 

Pe = [0; 0; 0] ; 

To = [15 ; 0 ; 0] ; 

*/.U.U. inertia tensor matrix 

lb = [ 10 0 0; 0 16.12 0; 0 0 7.97] 
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APPENDIX B: SIMULINK FILES 



The nine state nonlinear model of bluebird is contained in the SIMULINK 
file E0M.9.m and was trimmed at a flight condition of 

• flight speed equal to 73 fps 

• flight path angle equal to zero 

using the TRIM command. The resulting trim values for the state vector and 
input vector are: 



x 



73.3 

0 

-0.0023 

0 

0 

0 

0 

0 

0 



and u = 



0 

0 

0 

0.2858 



The LINMOD command was used to linearize the sixteen state nonlinear 
model of Bluebird (contained in the SIMULINK file EOM.16.m and described 
in Chapter III) about this trim point. The resulting linear system is contained 
in the MATLAB file Linear 16. mat. 

The rudder was removed as a control input (remove the second column 
of the B matrix) and the resulting linear model was used as a basis for the 
synthesis model contained in the SIMULINK file synthesis. m. This synthesis 
model was used to determine the LQR gain. The synthesis model, Q and R 



98 



weighting matrices, and resulting LQR gain is contained in the MATLAB file 

LQR.dat.mat. 

The full nonlinear simulation is contained in the SIMULINK file, plant l.m. 
The MATLAB .m file simdata loads the workspace with the appropriate vari- 
ables. The file simdata calls the .m file trajectory. m in order to generate the 
trajectory schedule. Any changes to the commanded trajectory or wind distur- 
bance schedule can be made in trajectory. m. 

A version of the nonlinear simulation that does not use the filter and 
sensor blocks is contained in the SIMULINK file plant2.m. It runs considerably 
quicker than plant l.m. 
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APPENDIX C: ©-IMPLEMENTATION 
PROOFS REFERENCED 

Identity. Let x € iZ 3 = const be given. Then 

^(ifl(A)x) = ~' b R( A)* x S-'(A). (C.l) 

and 

^(fi?(A)x) = x xf *(A)S-'(A). (C.2) 

Proof: To derive both equations we will need Poisson’s Law: 

£('bM A)) =" u,„ xi JH A), (C.3) 

and the following identity: 

a x b = —b x a (C.4) 

for any vectors a and b of compatible dimensions. Now, consider 

j t (' B R(A)x) = (1(’ b R( A))* +£ H(A)|* 

= B o;b Xg R(A)x = —gR(A)x x B Wfl, (C.5) 

using equation (C.3), (C.4) and x = const. 

Next, by the chain rule we get 

Tt^Ri A).) = ^(ilKA)»)|A 

= -L( B R(\)x)S(\) (C.6) 

Equation (C.l) now follows by comparing equations (C.5) and (C.6). 
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To obtain equation (C.2), consider 



j t (' B R(A)fR(A)) = R(A))f R(A) +' R(A)£(?R(A» = 0, (C.7) 

since gR(A)fR(A) = I VA. Now, using equations (C.3) and (C.7) we get: 

j ( (?R(A)) = -?R(A) b u, b x . (C.8) 

Finally, following the steps in the derivation of equation (C.l) we obtain: 

Jf(?*(A)*) = x xf R(A)S-'(A). 



Theorem ..1 Suppose that assumptions A\ through A3 hold. 

Al. Dim(x C 2 ) — dim(u) = dtm(y 2 ). 

A2. The matrix 

si — Ac i 2 ?c3 

~C cl C C 2 

has full rank at s = 0. 

A3. The matrix pair (A c ijC d) « observable. 

Then for each equilibrium point of Q in £ the following properties are 
observed: 

• the feedback systems %{Q , K )(r 0 ) and T(Qi 0 , Ki ) have the same closed 
loop eigenvalues; 



T,(g ,K )(r 0 )(s) = L(A 0 )T(Gi 0 , )C, )(s)L~ l ( A 0 ), 

Ao = [0 1 ) 1 ^ 
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Proof: In the proof we set the controller matrices D c i, D c 2 to zero. This 

does not change the results but considerably simplifies the algebra. Further- 



more, we will drop explicit dependence of the controller parameters on a. Let 

Pi b 



(x VQ ,x Po = ^ , Uo, ro € S) be given. Consider the feedback interconnection 

of the linear plant Qi 0 (A 0 ) and linear controller ACi . The state matrix F of this 
feedback system has the following form: 



F := 



Ax 


a 2 


B 2 C c1 


B 2 C c2 ■ 


I 


4 


0 


0 


B c iCi 


0 




Be 3 


c 2 


0 


0 


0 



(C.9) 



Next we linearize the feedback interconnection of the plant Q and the controller 
K . However, in order to that, first we must determine the values of the con- 
troller states x c i and X& along the trajectory r 0 € S. From equation (4.9) we 
obtain: 



^£ci 0 — A c ix c i 0 + B c i—yi 0 + B&L 1 (A 0 )eo 

d , 

^*c 2 0 = C c ix clo + C c2 L (A 0 )e 0 



Uo — X C 2 0 . 



Notice, since along r 0 : 



Co = 0, y i 0 = const , x c2o = u 0 = const 



we get 



^ ^clo — A c lXclo 

0 = C c \X c i 0 . 
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Now, using Assumption A3 we conclude that x c i 0 = 0. 

In order to compute the linearization of the feedback interconnection of 
Q and K, , we must first obtain the linearizations of equations (4.3) and (4.9) 
about the operating points (x^Xp^tto) € E and (x c i 0 = 0,x c2o = «o,yi 0 = 
0, eo = t/ 2 0 — ro = 0) respectively, determined by tq € E. The linearization of 
the plant Q is given by (4.4). The linearization of the controller K has the 
following form: 

£ci = Acifd + Bci 0i + B c3 B- 1 (0 2 — p) 

(c2 = C el U + C c iL-\d 2 — p) 

V = U (C.10) 



It is easy to verify that the state matrix M of %{Q , C )(r 0 ) is 



M := 



Ai Aj 0 B 

L A 4 0 0 

BeiCiAx B c iCiA2 + Bc 3B -1 Ad Bel Ci B 

0 CdL- 1 Cel 0 . 



Let 



Obviously, 



'7 0 0 0 
0 7-00 
0 0 7 0 
0 0 0 7 



P - 1 = 



'7 0 00' 

0 7" 1 0 0 
0 0 7 0* 

0 0 0 7 



Now using simple algebra it is easy to show that 



(C.ll) 



St = P-'MP 

Ai A 2 0 B 

7 A, 0 0 

BciCiAi B c i Ci A 2 + Bc 3 Ad BciCiB 

. 0 C c 2 Cel 0 
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A\ A2 
I A, 




0 

0 


B ‘ 
0 




B c i 


Ci 0 

0 I 


A\ A 2 
I At 


+ Bc 3[0 I] 


Ad 


B c , C\B 








0 C c 2 




Crt 


0 J 



The proof of the first part of the theorem now follows from Assumption A2 and 
an observation that the matrices F and M axe in the form of the matrices F 
and M. The proof of the second part of the theorem consists of the following 
steps: 



1. compute the transfer matrix of the linearization of the controller K (A) 
using equation (C.10) from controller inputs 0j, 62, p to controller output 
Vi 



2. apply a similarity transformation 

I 0 

0 L~ l ' 




to the linearization of the plant Q (= Qi (ro)) given by equation (4.4) and 
derive the transfer matrix from the control inputs of the linear plant 77 to 
the outputs 0\ and 62 using this new state-space realization; 



3 



compute the feedback interconnection of the transfer matrices obtained in 
steps 1 and 2 to get the final result. 



A simple computation shows that the transfer matrix from the controller inputs 
#1, 62, p to the controller output 77 is given by: 



«« = C c ,(sl - -(h(s) - fts)) 

s 

+ B c i <?j(s) + C C 2L~ l —(§2(3) — p(s)) 

s 





Pi 


■ «,(.) ' 


“ 


K,(,) 


h(s) 






L 




m 
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where Ki(s) is the transfer matrix for the controller £/ . 



Applying similarity transformation P\ to equation (4.4) and computing 
the transfer matrix from tj to 0\ and 0 2 results in: 



«t(s) j 



A\ A 2 


B ' 


I Ai 


0 


Ci 


0 


Di 


. 0 


L 


0 






(C.13) 



where the transfer matrix is given in packed matrix form. A simple observation 
shows that 



' <Ms) ' 

W J 



(#). 



where G r t,(s) is the transfer matrix for the plant Qi 0 . 



Now routine algebra shows that the transfer matrix from p to of the 
feedback interconnection of the transfer matrices in equations (C.12) and (C.13) 
is given by: 



T,(C X )(r„)« = L(A „)T(S h , K, )(«)£-•( A„). 
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